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NOMENCLATURE 


Symbol Definition 

A Reference area of particle 

Be Best number (Be = C^Re^) 

Cj Drag coefficient 

Steady-state drag coefficient 

Cda coefficient due to angle of attack 

Lift coefficient 

Steady-state lift coefficient 

Lift coefficient due to angle of attack 

Moment coefficient 

Steady-state moment coefficient 

Moment coefficient due to angle of attack 
ma 

D Particle diameter 

D Aerodynamic drag force 

E Total collection efficiency 

FRL Flight reference line 

g Accleration of gravity (9.8 m/s^) 

g^ Conversion factor (1 kg-m/newton-sec^) 

h Projected height of the body along the vertical coordinate 

line 

Moment of inertia of mass relative to the z axis 
L Aerodynamic lift force 

LWC Liquid water content 

Reference length of body 


Moment of aerodynamic forces acting on the particle 
Mass of water droplet (kg) 

Mass flow rate 

Number of particles of size i per unit volume 

Reynolds number based on the diameter of the particle 

Surface distance measured from the leading edge of each 
body; positive along the lower surface and negative along 
the upper surface 

Lower surface impingement limit 

Upper surface impingement limit 

Time 

Velocity of particle 

Velocity of particle relative to flow field 

Velocity of flow field 

x-component velocity of flow field 

y-component velocity of flow field 

Free-stream velocity 

x-coordinate of particle 

x-component velocity of particle 

x-component acceleration of particle 

Initial value of the horizontal coordinate of particle 

y-coordinate of particle 

y-component velocity of particle 

y-component acceleration of particle 

Initial value of the vertical coordinate of particle 

Upper tangent trajectory of the particle corresponding to S 


Symbol 


Definition 


Lower tangent trajectory of the particle corresponding to S 


Greek Symbols 

a Angle of attack 

3 Local collection efficiency 


Y 

6 

n 

e 

• 

e 

0 


Particle path angle 


Y = tan 



The angle between ^ and 

a 

Percentage liquid water contained in particles of size D^- 
Pitch angle of particle 
Angular velocity of particle 
Angular acceleration of particle 


1.0 INTRODUCTION 


# 


General aviation aircraft and helicopters exposed to an icing envi- 
ronment can accumulate ice resulting in a sharp increase in drag and 
reduction of maximum lift causing hazardous flight conditions. NASA 
Lewis Research Center (LeRC) is conducting a program to examine , with 
the aid of high-speed computer facilities, how the trajectories of 
particles contribute to the ice accumulation on airfoils and engine 
inlets. This study, as part of the NASA/LeRC research program, develops 
a computer program for the calculation of icing particle trajectories 
and impingement limits relative to axi symmetric bodies in the leeward- 
windward symmetry plane. 

The methodology employed in the current particle trajectory cal- 
culation is to integrate the governing equations of particle motion in a 
flow field computed by the Douglas axisymmetric potential flow program 
[1]. The three-degrees-of-freedom (horizontal, vertical, and pitch) 
motion of the particle is considered. The particle is assumed to be 
acted upon by aerodynamic lift and drag forces, gravitational forces, 
and, for nonspherical particles, aerodynamic moments. The particle 
momentum equation is integrated to determine the particle trajectory. 
Derivation of the governing equations and the method of their solution 
are described in Section 2.0. 

General features, as well as input/output instructions for the 
particle trajectory computer program, are described in Section 3.0. The 
details of the computer program are described in Section 4.0. Examples 
of the calculation of particle trajectories demonstrating application of 
the trajectory program to given axisymmetric inlet test cases are pre- 
sented in Section 5.0. For the examples presented, the particles are 
treated as spherical water droplets. In Section 6.0, limitations of the 
program relative to excessive computer time and recommendations in this 
regard are discussed. 
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2.0 METHODOLOGY 


The procedure for computing the particle trajectory around an 
axi symmetric body has been divided into the following steps; 

1. Compute the potential field around an axisymmetric body 
assuming the particles do not influence the flow field. 

2. Generate grids around the body to satisfy the refinement 
level and velocity error criteria which are input by the 
user. 

3. Determine the velocity flow field on the grids around 
the body with or without angle of attack for a given 
free-stream airspeed. 

4. Calculate the trajectories of droplets in the flow field 
determined in Step 3 using the Adams-Moul ton predictor- 
corrector method to integrate the equations of particle 
motion. 

5. Calculate the local collection efficiency for the body. 

The computational procedures used in Step 1 are described in 
Section 2.1, those used in Steps 2 and 3 are discussed in Section 2.2, 
and those used in Steps 4 and 5 are described in Section 2.3. 

2.1 Potential Flow About an Axisymmetric Body 

The Douglas axisymmetric potential flow program developed by Hess 
and Smith [2] is used in the present study for calculating the flow 
field about the body. This computer program uses a distribution of 
sources, sinks, and/or vortices along the body surface to calculate the 
potential flow field. The body surface is represented by an arbitary 
number of straight line (or curve) segments. In calculating the flow 
field, contributions from all the sources, sinks, and/or vortices are 
summed. The accuracy of this method was tested by comparing its pre- 
dicted velocities and surface pressure coefficients with both analytic 
solutions and experimental data [2]. Excellent agreement has been 
found. 
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The Douglas axisymmetric flow program consists of essentially three 
parts; a geometry-generating program called SCIRCL and an axisymmetric 
flow field computer program called EOD. The SCIRCL program generates 
the geometry input to the potential flow program (EOD) for a given 
specified analytical shape. Using the input from SCIRCL the EOD program 
calculates the flow field at any position in space. The EOD program is 
used to generate an input data tape containing the sources, sinks, and 
vortices which are then used to compute the velocity at each time step 
along the particle trajectory. 

Only limited details of the potential flow program are provided in 
this report since the major thrust of this study was to integrate the 
axisymmetric potential flow program into the program for computing 
particle trajectories and local collection efficiencies documented in [3]. 
Details of the particle trajectory program is reproduced in this section 
for completeness. The already-developed Douglas potential flow program 
is simply used to provide flow field input [1]. For complete details of 
the potential flow program, the user should consult References 1 and 2. 

2.2 Particle Equations of Motion 

In the present study the motion of a particle has been analyzed as 
a point mass particle which is acted on by the potential flow field but 
which itself does not influence the flow. The forces acting on the 
particle are considered to be those of lift, drag, and gravity. Pitch 
moments acting on the particle are also considered. Figure 2.1 shows 
the forces acting on the particle and the velocity vectors relative to 
the motion of the particle. The flight reference line (FRL) shown in 
Figure 2.1 is not significant for a spherical particle; however, for 
more arbitrarily shaped particles, e.g., a snow flake, the FRL must be 
defined relative to the lift, drag, and moment coefficient data avail- 
able for the given particle shape. The governing equations of the 
particle motion are derived for the more general case, i.e., arbitrary 
particle shapes. In turn, the computer program developed in this study 
and described in Sections 3.0 and 4.0 provides the option for general- 
shaped particles. However, the test cases given in Section 5.0 are for 

spherical particles only. A valid drag law for spherical particles is 
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Figure 2.1 Diagram of the velocity vectors and forces acting on a point 
mass particle. 

built into the computer program. The user must input lift, drag, and 
moment coefficient data for general -shaped particles. 

The equations of motion of the particle derived from a force 
balance on a point mass particle as shown in Figure 2.1 are: 

mx = -D cos Y-Lsiny (2.1) 

my = -D sin y + L cos y - nig (2.2) 

where 

1 y - W, 

Y = tan” ^ (2.3) 

The flow field velocity components in the longitudinal and radial direc- 
tions, i.e., W and W , respectively, are obtained from the potential 
X y 
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flow program described in Section 2.1. The aerodynamic drag and lift 
forces are defined as: 


D = 


Pa^a^ 

^d "2^^ 


L = 


C ° ° A 
2g^ ^ 


where A is a characteristic area of the particle, p. is the density of 

a 

air at the position of the particle, and V, is the particle velocity 

a 

relative to the flow field and is defined as: 


Vg = /(X - W^)2 + (y - W^)2 

(2.4) 

For arbitarily shaped particles, expressions for the drag 

coefficient. 

C^, and lift coefficient, must be provided by the user 

often approximated with: 

. They are 

^d " ^do ^ ^da“ 

(2.5) 

'4 ' '^40 * '^4«“ 

(2.6) 


where a is the angle between the FRL and the velocity vector (Figure 

a 

2 . 1 ). 


Equations 2.1 and 2.2 can be solved given values of the coefficients 
of Equations 2.5 and 2.6, i.e., ^£a’ angle a is 

computed from the expression: 

a = 0 - Y (2.7) 

where the angle e, generally called pitch angle in aerodynamics, is 
governed by: 

0 = ^ (2.8) 

^zz 

where is the moment of inertia of mass relative to the z axis. The 
moment of aerodynamic forces acting on the particle is: 


where is approximately: 

C = C + C *a 
m mo ma 


and are constants to be provided by the user depending upon the 
shape of the particle, and is a reference length. 

For a spherical particle with zero angular velocity, the lift force 
is always zero for potential flow. The governing equations are signi- 
ficantly simplified in this case: 

mx = -D cos Y (2.10) 

my = -D sin Y - mg (2.11) 

where 



Pa^'a' 

2g. 


( 2 . 12 ) 


In the present study is the steady-state drag coefficient of a 
sphere as a function of Reynolds number based on particle diameter, Re = 
V^D/v, as shown in Figure 2.2. The diameter of the spherical particle, 

G 

D, and the kinematic viscosity of the air, v, are assumed constant along 
the particle's trajectory in the present study. The drag law utilized 
was provided by NASA/LeRC (Figure 2.2 and Table 2.1) [5]. 

The governing Equations 2.1, 2.2, and 2.8 together with the follow- 
ing definitions form a complete set of equations to describe the motion 
of particles in the flow fields. 


dx _ • 
dt ^ 

(2.13) 

II 

(2.14) 

de _ A 

dt ■ ® 

(2.15) 


Integration of the Equations 2.1, 2.2, 2.8, 2.13, 2.14, and 2.15 results 




TABLE 2.1 Polynomial Coefficients Relating Best Number (Be) to Reynolds 
Number for Spherical Particles [5]. 


n 

Be = 1 

j=o 

a.Re^ 

u 


Reynolds Number Range 

j 


0.05 < Re < 3 

0 

0.0 


1 

24.167 


2 

3.254 


3 

-0.23564 

3 < Re 5 330 

0 

-28.339 


1 

38.969 


2 

0.73204 


3 

-0.00056084 

330 < Re 

0 

0.0 


1 

93.462 


2 

0.37576 


where 


Re = 

V 


Be 

D 

-> 

W 

V 


Best number (C^Re^) 
particle diameter 
inertial velocity of the air 
(W = Wj + Wyj) 

inertial velocity of the particle 
(V = xi + yj) 
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in solutions for x, y, e, x, y, e of the particle at time t. For the 
sphere particles for which the lifting force is omitted, only Equations 
2.10, 2.11, 2.13, and 2.14 are integrated. The integration method is 
described briefly in the next section. 

2.3 Particle Trajectory Solution Algorithm 

The method utilized for integrating the governing equations of 
particle motion is the Adams-Moulton predictor-corrector method [6]. 

The solution is obtained if the summation of the difference between the 
particle velocities in the x and y directions obtained by the predictor 
and corrector, respectively, and divided by the value of the solution 
at the last time step is less than a specified number e: 


(X - X ) 

n, predictor n, corrector^ 


(9 - Y ) 

' n, predictor n, corrector 

*n-l 

^-1 



The value of e is specified by the user. 


2.4 Computation of Collection Efficiency 


Particle trajectories calculated as described in the previous 
sections are used to establish the relations between the particle's 
initial position (-«,yQ) and the position where it impinges on the body 
surface, s. s is the length along the body surface measured from the 
leading edge on the body to the point of particle impingement. The 
value of s is defined as positive on the lower surface and negative on 
the upper surface, y^ is the initial value of the vertical coordinate 
from which the particle is released (see Figure 2.3). The local collec- 
tion efficiency, 6, is calculated as a function of the distance along 
the body surface by differentiating y^ = y^is) with respect to s [7]; 



(2.16) 


The minus sign is introduced so that 6 is positive, which is consistent 
with the definition given in Reference 4. 


Upper surface impingement limit 



Figure 2.3 Illustration of impingement terminology and water droplet trajectory in an 
airfoil flow field [ 7 ]. 




The overall collection efficiency, E, is defined as: 



where and y^^ are the upper and lower tangent trajectories of the 
particle relative to the body surface and h is the projected height of 
the body along the vertical coordinate line. 

2, 5 General Computational Procedure 

The general computational procedure to determine the local collec- 
tion coefficient of a body is carried out as follows. First, the 
initial conditions for the differential equations governing the particle 
motion are determined either automatically by the computer program or 
input manually by the user. These conditions call for specification of 
an initial particle position x^,y^ and an initial particle velocity 
^o’-^o‘ computer program described in the following sections auto- 

matically determines the initial conditions if desired. The initial 
upstream x-coordinate, x^, is assigned the value of x at which the 
difference in the free-stream velocity W and the local velocity W is 
less than some small value e. The value of e is specified by the user. 

To determine the initial vertical coordinate y^, the computer 
automatically searches for the upper and lower limits of the y- 
coordinate, y^^ and y^^, respectively (see Figure 2.3). Any particles 
released within this region will strike the body. Any particles released 
outside this- region will miss the body and are of no interest to the 
computation of collection efficiency. The range of vertical position 
you to y^^ is then divided into a number of increments prescribed by the 
user. The trajectory of particles leaving each of these vertical posi- 
tions is calculated and the impingement position of the particle on the 
body surface, s, is recorded. This collection of (Yq.s} values plus 
those generated during the computer search for y^^ and y^^ are used to 
express s as a function of the particle's initial vertical coordinate 
y^. The value of 6 = -dy^/ds is then computed by a linear approximation 
or by curve fitting the total collection of data points (yQ>s} to a 


polynomial curve fit. The degree of the polynomial is specified by the 
user. The current program allows the user to curve fit the entire set 
of datum points to one curve or to fit the curve in segments using a 

prescribed number of points on either side of the specified position 
(see Figure 2.4). This latter procedure is similar to segment-averaging 
or segment-curve fitting of the entire curve. 

The initial velocity of the particle is prescribed to be equal to 
the value of the flow field at x^,y^, i.e., x = and y =, 

'^y(Xo»yo) otherwise specified by the user. 

The following sections describe in detail the computer program and 
necessary user's information to compute particle trajectories and local 
collection efficiency for two-dimensional airfoils and inlets. 


Initial y Position of Particles, Initial y Position of particles, y /£ 

• . .O.C I I iO|C 


x-x- 


Input data 
Curve fit 


■0.13 1 1 1 1 1 1 1 1 

-0.04 -0.02 0.0 0.02 0.04 0.06 0.08 0.1 ( 

Distance along Body Surface, s/t^ 

(a) Total data set curve fitting; 6 degree polynomial with all input 
data. (Joukowski airfoil 0015, a = 4“, with walls). 




X : Input data 

: Curve fit 


Distance along Body Surface, s/t 


(b) Segment curve fitting; 3 degree polynomial with sequential sets 
of input data. (Joukowski airfoil 0015, a = 4^^, with walls). 

Figure 2.4 Illustrates total data set and segment curve fitting techniques 




3.0 GENERAL DESCRIPTION OF COMPUTER PROGRAMS 


The purpose of this section is to describe the computer program in 
sufficient detail so that it can be run successfully by the user. 

Section 3.1 describes some general features of the program which will 
better enable the user to follow the data input instructions given in 
Section 3.2. Instructions for the geometry generation program and for 
the axisymmetric potential flow program are also given in Section 3.2. 
These, however, are simply reproduced from Reference 1 without appre- 
ciable discussion. Input instructions for the grid generation computer 
code are also given in Section 3.2. The original references should be 
consulted if additional information is required. 

3.1 General Features of the Program 

3.1.1 Types of Flow 

Axisymmetric potential flow over arbitrarily shaped bodies is 
considered. The water droplet in the flow field is treated as a solid 
sphere although the option for a nonspherical particle is provided in 
the computer program. For nonspherical particles the user must provide 
expressions for the coefficients of aerodynamic lift, drag, and moment. 

3.1.2 Surface Distance Computation 

The surface distance along the body is computed by summing elements. 
As, determined from a linear approximation, as = /(ax)2 + (Ay)2. 

Surface distance is measured from the leading edge of the body with 
positive values defined on the lower surface and negative value defined 
on the upper surface. Figure 2.3. 

3.1.3 Initial Longitudinal Coordinate, x^, of the Particle Trajectory 

The initial longidinal coordinate position x^ at which the particle 
trajectory calculations begin (Figure 2.3) is automatically determined 


by the computer program. The value of is selected by testing the 
maximum difference between the locally computed value of W and the 
free-stream velocity at successively farther distance upstream. The 
value of X for which the inequality 

< 0.001 (3.1) 

max 

is satisfied is designated as x^. Equation 3.1 is tested over a speci- 
fied range YYLO < < YYUP, where YYLO and YYUP are initial values 

input by the user. They represent the expected maximum and minimum 
range of the upstream y-coordinate. ''s the reference length which is 
normally the chord length for an airfoil or the mouth diameter for an 
inlet. This procedure for selecting the initial position of the particle 
trajectory has been developed so that computer time may be conserved by 
starting the particle as near the leading edge as variations in flow 
field velocity will allow. 

3.1.4 Impingement Position of the Particle on the Body 

A coordinate transform technique is utilized in determining the 
position at which a particle strikes the body. The transform technique 
is illustrated in Figure 3.1. The equations governing the coordinate 
transform are: 



Xt = X 

Y _ V - VREF 

’t YREF - YLO 

where (see Figure 3.1a) 

YREF = if Y > 0.0 


(3.2) 


YREF = Y^-ab-CC- Y < 0.0 

where X and Y are the x- and y-coordinates normalized by the characteris- 
tic length, 2,^, respectively. If the particle moves across the coordinate 
line Y^ = 0 in the transformed plane and the position of the particle 
is greater than zero, the particle is recorded as having crossed the 
surface of the body. An iteration procedure, described in the following 


X = XIN 



Note : XREAR is initial input data, YUP and YLO are automatically 

calculated in the program. 

(a) Physical plane 



Figure 3.1 


Coordinate transform. 
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subsections, is then carried out to determine the exact surface location 
of impingement. 

3.1.5 Computation of Surface Impingement Location 

The method of computing the surface location of particle impinge- 
ment is described in this section. First a general case and then a 
special case (see Figure 3.2) are considered. 

3. 1.5.1 General Case . During a time step At, consider the particle 
to cross the body along the line ON (see Figure 3.2a). Let the coordi- 
nates of the particle at time t be (X0,YP0) and the coordinates at t = 
t + At be (XN,YPN). The reference coordinates on the body surface 
A'ABCC or A'AB'CC are X0,YR0 and XN,YNR, respectively. The analytical 
function of the line ON" is: 


Y - YPO YPN - YPO 

X - XO XN - XO 


The function of the line O'N', which joins the two nodal points describ- 
ing the body surface, is: 

Y - YRO YRN - YRO /o c 

X - XO ■ XN - XO ^ 


The coordinates (XI, YI) of the intersection of line ON and O'N' found by 
simultaneous solution of Equations 3.4 and 3.5 are: 


YT (XN - X0)(YP0 - YRO) , y. 

“ (YRN - YRO - YPN + YPO) 

YT - (YPN - YPO) (YPO - YRO) , yp. 
^ “ (YRN - YRO - YPN + YPO) 


(3.6) 


The reference coordinates on the body surface (XI, YR) corresponding to 
the position (XI, YI) is then found. If |YI - YR| < 10 the point 
(XI, YR) is regarded as the particle impingement position on the body 
surface. If |YI - YR| > 10”^ the point (XI, YI) is redefined as (XN,YPN) 
and (XI,YR) as (XN,YRN), respectively. The procedure is repeated until 

_C 

the requirement |YI - YR| < 10 is satisfied. 


3. 1.5. 2 Special Case . Consider the particle crossing the body 
surface for the special trajectory shown in Figure 3.2b. The function 
of line OFT is: 


. 

‘ F-r 


Particle — 
Trajectory 
(Special Case) 


\ B' 

Particle 
Trajectory 
(General Case) 


Particle Trajectory 


^ (X0,YF 


(XI, YR) 


Body Surface 


r Boay 


(Xri,YRN) 


(XI, YI) 1 (XN,YPN) 


i/(X0,YR0) 

/ 

Time Step (N-1) Time Step (N) 
(a) General case 


Body Surface 


J'^(XN,YRN) 

Particle Trajectory-i , . yy ^ 

A (XF,YF) [ 

° (X0,YP0) 

Time Step (N-1) Time Step (N) 

(b) Special case 

Figure 3.2 Illustration of the method of determining the position 
of particle impingement on the body surface. 
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Y - YPO YPN - YPO 
X - XO " XN - XO 


(3.7) 






and of 1 ine AN ' is 


Y - YF _ YRN - YF 
X - XF “ XN - XF 


(3.8) 


The coordinates (XI, YI) of the point of intersection of line ON and line 
AN' are found as follows: 


Let 

A = (XN - XF)(YPN - YPO) - (YRN - YF)(XN - XO) 

B = YPN - YPO 
C = YRN - YF 
D = XN - XF 
E = XN - XO 

then 

XI = [B-D-XO - E-C-XF + (YF - YPO)-E*D]/A 

YI = [D-B-YF - E-C-YPO + (XO - XF)-B-C]/A (3.9) 

The same procedure and criteria as for the general case is used to 
determine the position of particle impact on the body surface. 

3.1.6 Upstream y-Coordinate Limits 

The method by which the computer selects the upper and lower trajec- 
tory limits y^^ and y^^ is described in this section. A number of 
options are available to the user depending on the setting of the flags 
LRANG and LIM. If LRAN6=1 and LIM=1 , the program searches automatically 
for the upper and lower limits of the radial (vertical) coordinates of 
the initial particle position y^^ and y^^ (see Figure 2.3). 

The search procedure consists of the computer program initially 
seeking the range within which y^^^ and y^^ lie. This is achieved by 
computing the particle trajectory from the initial position (X0,YMAX) 
where YMAX and YMIN are the user's initial guessed values of y^^ and 
y^^, respectively, and XO = The particle passes under or 

hits the body, YMAX and YMIN are redefined as YMAX^^^^^ = YMAX + aY and 
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VMINn^I^ = YMAX where aY = YMAX - YMIN. The procedure is then repeated 
until the particles pass over the body. The current values of YMAX and 
YMIN then specify the range within which more precise values of y^^ and 
are sought. If on the first trajectory calculation the particle had 
passed over the body then the trajectory from (X0,YMIN) is computed. 

If the particle again passes over the body then the above procedure is 
reversed (i.e., YMAX,^^^ = YMIN and YMIN^^^^ = YMIN - AY). 

Once this order of magnitude range of y^^ and y^^ is determined, 
more precise values of these limits are computed as follows. A particle 
trajectory from the position Y' = (YMAX + YMIN)/2 is computed. If the 
particle passes under or hits the body then the next trajectory is com- 
puted from = (Y' + YMAX)/2. Alternatively, if it passes under 

the body then = (Y' + YMIN)/2. Successive halving of the range 

YMAX to YMIN in this manner continues until convergence is achieved. 
Convergence is assumed when the difference of the y^ coordinate between 
two trajectories, for which one impinges on the body and one misses the 
body, is less than a small value specified by the user. In the program 
the small number is designated as YLIM. For the test cases given in 
Section 5.0, YLIM = 10“ m was used. Determination of the values of y^^ 
and y^^ also provides values of the upper surface impingement limit, Sy , 
and the lower limit, Sy , respectively. 

If the control flags are set such that LRANG=0 and LIM=1 , the 
computer program will search for y^^ and y^^ within the range YMAX and 
YMIN input by the user. If a poor guess was made and y^^ and y^^ are 
not in this range, then the program is terminated. If the user desires 
to compute only one particle trajectory, the flags should be set to 
LRANG=0 and LIM=0. 

3.1.7 Calculation of the Local Collection Efficiency 

With the values of y^^ and y^^ determined, the computer divides 
the range into NPL segments (i.e., (y^^ - yQj^)/MPL) which represents the 
number of particles to be computed. Particles are then released from 
each of the NPL locations. The surface locations of impingement of each 
particle is then recorded and a collection of coordinates IYqjS} are 
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stored in a file. These values are used to construct the functional 
relationship between and s from which the local collection efficiency 
is calculated. 

Two methods are used to calculate the local collection efficiency, 
6. The first method is linear approximation; 


The second method utilizes a polynomial curve fit. The computed values 
of surface impingement, s, are determined as a function of y^, i.e., 
s = sCy^), by a polynomial curve fit. The value of 6 is computed by 
taking the derivative of the polynomial function. The number of coeffi- 
cent of the polynomial function is input by the user as NCOEF (the order 
of the polynomial is then NCOEF-1). The total number of points (YqjS} 
can be curve fit or the function can be curve fit in segments similar to 
a running average. For segment curve fitting the variable NS=0 and the 
number of sequential coordinate points used in a segment is specified by 
the user with the variable NPTS. If NS=1 , all points are used. For 
multi -size distribution cases (NSI>1), only the all-points polynomial 
curve fit option is available. 


3.2 Program Input Instructions 

A general flowchart for the computer program is shown in Figure 
3.3. The general procedure for running the program consists of first 
manually constructing a geometry input data tape. Tape 05. The geometry 
is specified in terms of a number of coordinate positions on the surface. 
The origin for the Cartesian coordinate system used to specify the 
geometry can be selected arbitrarily. The geometry generation computer 
program, SCIRCL, is then run to create Tape 17 which is the input file 
for the axisymmetric potential flow program. Section 3.2.1 describes 
the input tape or card deck structure for the geometry generation 
program. The format of the output tape. Tape 17, created by SCIRCL is 
the input to the axisymmetric potential flow program and is described 
in Section 3.2.2. 


21 


Tape 05 
Geometry 
(Input Data) 


Tape 02 
Particle 
(input data) 



3.3 General flowchart. 
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Now using Tape 17 as input, which must be renamed Tape 05, the 
axisymmetric flow code is run and creates Tape 21. Tape 21 is directly 
provided to the particle trajectory computer program. Additionally, 
Tape 02, which must be manually created as described in Section 3.2.3, 
is required to run the particle trajectory program. 

Running the particle trajectory program creates the output tapes 
listed below: 

Tape 01: Data stored on Tape 01 is later written to Tape 03. 

Tape 01 = Tape 03 for unsymmetric flow cases. 

Tape 03: Data stored on Tape 03 is used for the calculation 

of collection efficiency. 

Tape 04: Tape 04 is created if LOPT^^O and contains coordinates 

of the particle trajectories as well as values of 
the wind speed component as a function of time. 

Tape 08: Data sotred on Tape 08 is used for trajectory plots. 

The data sotred are the ,y/ 2 ,(, coordinates of the 
particle XP,YP for each time step and the surface 
impingement point SW. SW is initially set at 
88.8888. When SW;^88.8888, the particle trajectory 
is terminated. 

Tape 09: Tape 09 contains the local collection efficiency, B, 

as a function of surface position, s/i^. The program 
automatically plots 6 versus s/ic at the NASA/LeRC 
facility. 

The input instructions for running the computer programs are 
presented in the following sections. The instructions are given in 
terms of an input card deck structure which compares directly with data 
tapes or files. 


3.2.1 Geometry Generation Program Input Instructions 


The input instructions for the geometry generation program are 
taken directly from Reference 1. These instructions are given for 
continuity of this report. It is assumed, however, that the user is 
familiar with the details of the geometry generation program. The input 
deck for the geometry generation program has the following card structure 

Card Description 

1 

2 

3 

4 

5 (only if flag J > 0) 

6 

7 

8 1 


Number of '8' cards = NRAKES 


8 j 

9 

10 ■ 
11 
12 


10 

11 

12 . 


Number of '9' cards = ANBDYS 


Number of '10-11-12' groups for each '9' card = ANSEG 

*If ENREED = 99 on card 10, use 11a and 12a instead of 11 and 12 

*If ANSEG = 0 and TYPBDY ^ 0 on card 9, skip 10 and substitute 
11a and 12a for 11 and 12 
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Card 

No. 

Column 

Code 

Routine/ 

Format 

Explanation 

1 

1-36 

TITLE 

Main/ 

9A4 

Description of case 

2 

1-10 

XX 

Main/ 

F10.2 

Length, in plot inches, of 
x-axis required 


11-20 

XMIN 

tl 

Value, in data inches, of 
far left x point 


21-30 

EXEP 

II 

Data in per plot inch 
along x-axis 


31-40 

YY 

II 

Length, in plot inches, of 
y-axis 


41-50 

YMIN 

II 

Value, in data inches, of 
bottom y point 


51-60 

ORD 

II 

Data inch per plot inch 
along y-axis (usually 
equal to EXEP) 


61-70 

ELREF 

II 

The X values in area output 
data are nondimensional ized 
by ELREF. Default value is 1 


61-80 

AREF 

II 

The areas in area output 
data are nondimensional ized 
by AREF. Default value is 1 

For 

information to 

be passed on 

to EOD: 


3 

1 

IGEOMF 

Main/Il 

=1 , Use flat elements 
=0, Use curved elements 


2 

ISIGF 

II 

=2, Use constant source 
densities 

=1, Use linear source densities 
=0, Use parabolic source 
densities 


3 

ICURVN 

II 

=1, Read in curvature values 
=0, EOD will compute 
curvatures 


4 

NONEWF 

II 

=1, Use old velocity formula 
=0, New formula 


5-14 

ALPHER 

Main/ 

F10.2 

25 

If only one body is input, 
ALPHER is the angle of 
attack (used by EOD) 



Col umn 

Code 

Routine/ 

Format 

Explanation 

15 

IVORT 

Ma 1 n/ I 1 

= 1, Perform axi symmetric , 
closed-duct solution 
=0, Perform strip vortex 
solution 

16 

I PAR 

II 

Element geometry flag used 
by EOD 

=1, Parabolic elements 
=0, Linear elements 

17 

IFST 

II 

First-order terms flag 
=3, Both first-order terms 
=2, Curvature terms 
=1, First derivative terms 
=0, No first-order terms 

18 

ISND 

ll 

Second-order terms flag 
=3, Both second-order terms 
=2, Curvature squared terms 
=1, Second derivative terms 

19-20 

IFLLL 

Main/12 

=1, A combination solution 
will be calculated by 
EOD 

=0, No combination solution 
will be calculated by 
EOD 

1-8 

IDENT 

Main/A8 

Eight-character tag for 
case I.D. 

9-12 

PRC3 

Main/A4 

Title "EOD" 

17-20 

N06 

Main/14 

=0 

21 

LPNCHO 

Main/Il 

Flag A 

= 1 , Do not save output for 
EOD on Unit 17 

22 

IPLOTA 

Main/Il 

Flag B. Plot area against 
X position (see Reference 1) 

23-24 

IPLOTC 

Main/ 12 

Flag C 

=+l , Plot curvature versus S 
=-l , Plot curvature versus X 

25-26 

IREAD 

Main/12 

Flag D 

=0 (obselete) 


Card 

No. 

Col umn 

Code 

Routine/ 

Format 

Explanation 

4 

37 

lAB 

Main/Il 

Flag J. Redo geometry from 
point (see Reference 1) 


47 

IREDON(l) 

II 

Flag E. Redo entire geometry 
via direct interpolation 
(see Reference 1) 


48 

I REDON (2) 

II 

Flag F. LPNCHO for any redo 


49 

IRED0N(3) 

II 

Flag fi. IPLOTA for any redo 


50-51 

IRED0N(4) 

Mai n/I2 

Flag H. IPLOTC for any redo 


52-53 

IRED0N(5) 

II 

Flag I. TREAD for any redo 

All flag 

are on 

when equal to 1 

unless otherwise noted. (Either E or J or 

nei ther 

can be 

on but not both.' 

) Skip card 

5 if J = 0. 

5 

1-12 

XAA 

Main/ 

F12.5 

X position of starting 
point for partial redo 


13-24 

YAA 

II 

y position of starting 
point for partial redo 


25-36 

XBB 

II 

X position of ending point 
for partial redo 


37-48 

YBB 

II 

y position of ending point 
for partial redo 

6 

1-10 

ANBDYS 

Main/ 

F10.2 

Number of bodies 


11-20 

DELS 

II 

Spacing between points in 
region of interest 


21-30 

DELSMX 

II 

Maximum spacing far from 
region of interest 


31-40 

XRI 

II 

Axial distance at which 
surface distance equals 
zero 

7 

1-4 

NRAKE 

Main/14 

Number of axial rake locations 

8 

1-8 

XRAK 

Main/ 

F8.5 

Axial location of rake 


9-16 

YLO 

II • 

27 

y value of first point 
(lowest point) on rake at 
XRAK 


27 



Card 

No. 

Col umn 

Code 

Routine/ 

Format 

Explanation 

• 

8 

17-24 

YHI 

Main/ 

F8.5 

y val ue of last point 
(highest point) on rake at 
XRAK 

• 


25-27 

NY 

Main/ 1 3 

Number of points in rake at 
XRAK: Restriction INY<200. 

Rake points are equally 
spaced, AY, between YHI and 
YLO where 

• 





_ YHI - YLO 

- -(ny - 1) 



28-35 

XTRAN 

Main/ 

F8.0 

Value of axial translation 
of rake 

• 


36-43 

YTRAN 

11 

Value of vertical translation 
of rake 



44-51 

XSCALE 

II 

Value of axial scaling of 
rake 

• 


52-59 

YSCALE 

II 

Value of vertical scaling 
of rake 


9 1-10 TYPBDY Main/ 

FIO.O 


11-20 ASNEG 


Body number. However, if 
there is symmetry, then any 
body can be input as a 
mirror image of any other 
body. That can be accom- 
plished by setting TYPBDY= 
-M.N where M is the number 
of the body to be created 
and N is the number of the 
body to be copied. ANSEG 
is set to the Y value of 
the line about which body N 
is to be mirrored. No other 
input is required for this 
body except for ANLF. 

=Number of segments for the 
particular body, except as 
stated in TYPBDY 
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Card Routine/ 


No. 

Col umn 

Code 

Format 

Expl anation 

9 

21-30 

DELNEW 

Main/ 

= -l , Delta S spacing is set 




FIO.O 

to original value of 


DELS 

=0, Delta S is set to value 
of DELS from previous 
body 

=+number, Delta S is set to 
value of input 
DELNEW 


31-40 

ANLF 

II 

=0, Body is a lifting body, 
i.e., in EOD a vorticity 
solution about this body 
will be calculated 
=1, Body is a nonlifting 

body, i.e., no vorticity 
solution will be 
calculated 

Note: All lifting 

bodies must 

be input 

prior to any nonlifting bodies. 

41-48 

XTRAN 

• 

Main/ 

F8.0 

Value of axial translation 
of this body 

49-56 

YTRAN 

II 

Value of vertical trans- 
lation of this body 

57-64 

XSCALE 

II 

Axial scaling factor 

65-72 

YSCALE 

II 

Vertical scaling factor 

73-80 

XTMAX 

II 

Maximum value of x for 


which scaling is to be 
applied 

(Code indicating type of curve to be fitted through given points.) 

10 1-10 ENREED Main/ =0., for bisuperel 1 ipses [1]. 

F10.2 =1,000. Same as =0 but with 

finer point spacing near one 
end of segment (two such 
segments required). Usually 
used to give finer spacing 
at the highlight. The super- 
ellipse going into the high- 
light and the one coming out 
should have this flag. For 
bisuperel 1 ipses where the 
'1,000' option is to be 
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i 


Routine/ 

Format 


Explanation 


Main/ 

F10.2 


used, the rate at which the 
point spacing, ds , changes 
near one end dS . = dS - - 

(Rate) (dS^ _i ) can be ' ~ 
specified on input. 

The rate (program name 
= PACE) is entered as the 
fractional part of ENREED 
for each segment. For 
example, if ENREED were input 
as 1,000.06, the spacing for 
consecutive points would be 
evaluated as follows: 

DS. = DS._.| - (0.06)DS._.| 

if segment is to go from 
large to small spacing or; 

DS. = DS._.| + 1 .5(0.06)DS._.| 

if segment is to go from 
small to large spacing. 

If PACE is entered as 
zero (i.e., ENREED = 1,000.), 
the default value, 0.05, is 
used. (PACE < 0.133) 

*The first '1,000' super- 
ellipse ON A BODY reduces 
the point spacing as far as 
possible, down to a limit 
of 2 percent of the ds value 
at the beginning of the 
segment . 

*A11 subsequent '1,000' 
superellipses input will 
increase ds as far as 
possible up to the input 
value of OELS. 

*Any number or types of 
segments may be input 
between the first and 
subsequent '1,000' bisuper- 
ellipses, with the exception 
of a normal bisuperell ipse 
(ENREED=0.). 


Explanation 


= 1 . , Is a straight line, 
input 2 coordinates 
(XIN(l), YIN(l), XIN{2), 
YIN(2)). 

The first and last straight 
lines on bodies 2 and 3 and 
the last straight line on 
body 1 will automatically 
have their spacing increased 
from approximately DELS near 
the region of interest to 
approximately DELSMX away 
from the region of interest. 
To get this type of spacing 
in the first straight line 
of body 1 , ENREED must be 
specified as 10. 

=10., Special straight line 
used for initial straight 
line on lower shroud. The 
straight line starts with 
large spacing (DEMSMX) and 
ends with small spacing 
(DELS). 

= -l . , Fits a lemniscate 
between a straight line 
and a point. Input is 
three coordinates. 

=-3., Fits a cubic between 
two straight lines. Input 
4 coordinates. 

=-4.0, Generates a segment 
which is a mirrored image 
of all the points from 
(XIN(l), YIN(D) to (XIN(2), 
YIN(2)) about the line Y = 
YIN(3). See cards 11 and 
12 for XIN and YIN formats. 

=99., For direct interpola- 
tion option over one segment 
(see input instructions for 
card 12). 


Card 

No. 


10 


11 


12 


Column 

Code 

Routine/ 

Format 

Expl anation 

11-20 

REEDEN(l) 

Main/ 

F10.2 

Input exponent of x-term for 
bisuperellipse equation. 
Blank for all other segment 
types [1]. 

21-30 

REEDEN(2) 

II 

Input exponent of y-term for 
bisuperellipse [1]. 

1-12 

XIN(l) 

Main/ 

6F12.5 

x-coordinate for specified 
points 

13-24 

XIN(2) 



25-36 

XIN(3) 



37-48 

XIN(4) 



49-60 

XIN(5) 



61-72 

XIN(6) 



1-12 

YIN(l) 

Mai n/ 
6F12.5 

y-coordinate for specified 
points 

13-24 

YIN(2) 

If 


25-36 

YIN(3) 

II 


37-48 

YIN(4) 

II 


49-60 

YIN(5) 

II 


61-72 

YIN(6) 

II 



Note: If ENREED=99, input the following cards instead of cards 11 and 12. 


11a 


Z(l) Name list/ z is a complex array con- 

$B0DYIN/ taining the x value (in the 
real part) and y value 
(imaginary part) of each 
given point along the 
segment. The name list 
will normally be longer 
than one card. The program 
will use the input points 
with finer point spacing 
near regions of high 
curvature. 
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Card 

No. 

Column 

Code 

Routine/ 

Format 

Explanation 

12a 


DONE 

Name list/ 
$AUXIN/ 

A logical variable which 
should be input as: 
=.TRUE. 



BYPASS 


=.TRUE. if no refinement of 
points is required 


Note: If ANSEG=0 and TYPBDY^O, skip card 10 and substitute 11a for 11 

and 12a for 12. 



3.2.2 Axisvmmetric Potential Flow Program Input Instructions 


The input instructions for the axisymmetric potential flow program 
are taken from References 1 and 2 with some modifications. These instruc- 
tions are given for continuity. It is assumed, however, that the user is 
familiar with the details of the axisymmetric potential flow program. 

The input format given is the format of Tape 17 (renamed Tape 05 when 
input to EOD). The card or tape structure is as follows: 


Card 

Type Description 

1 Body title and case number 

2 Control flag card 

3 Chord/Mach number card 
4l 

Body control card 

6,6'] Input body element coordinates . 

6a [ IF0RMT=0 input card 6 

6b j IF0RMT=1 input card 6a 

IF0RMT=2 input card 6b 

7 Input curvature values (needed only if ICURVNj^O) 
Repeat cards 4-7 (NB+FLG05) times. 

8 Rake number card (needed only if IRAKE^O) 

9 Rake definition card (needed only if IRAKE^^O) 


Subroutine 

PARTI 

PARTI 

PARTI 

BASICl 

BASIC! 


BASICl 

BASICl 
I RAKES 


Repeat card 9 NN times. 

10 Nonuniform flow control flag (needed only if NNU^O) BASIC2 

11 Nonuniform flow velocities (needed only if NNU^O) BASIC2 

Repeat cards 10 and 11 NNU times. 


i 
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Card 

No. 


1 


2 


Column 

Code 

Routine/ 

Format 


Explanation 

1-60 

HEDR 

PARTI/ 

15A4 

Body description. 

66-69 

CASE 

PARTI/ 

A4 

Case number. 

1 

NB 

PARTI /II 

Number of bodies (I<NB<9). 

2 

NNU 

It 

Number of nonuniform flows (0<NNU<5). 

3 

FLG03 

It 

Ax i symmetric flow flag. 

4 

FLG04 

II 

Crossflow flag. 

5 

FLG05 

11 

Off-body point input flag. 

6 

FLG06 

II 

Basic data only flag. 

7 

FLG06 

II 

Ellipse generator flag (consistent 
with card 5), 

a 

FLG08 

II 

\ 


9 

FLG09 

II 


• Blank 

10 

FLGIO 

II 

J 


11 

FLGll 

II 

Perturbation velocities only. 

12 

FLG12* 

U 

Solve potential matrix. 

13 

FLG13 

II 

Blank. 

14 

FLG14 

II 

Prescribed tangential velocity. 

15 

FLG15 

II 

Strip-ring vorticity flag. 

16 

FLG16 

It 

Omit axisymmetric uniform flow 
solution. 

17 

FLG17 

It 

Omit crossflow uniform flow solution. 

18 

FLG18 

II 

Surface vorticity (instead of 
sources) for the final bodies. 

19 

FLG19 

11 

Prescribed values of the surface 
vortex strengths for the final bodies 
will be input. 


*Avai1able if and only if N0NEWF=1 , ISIGF=1 , and IGE0MF=1 (see card 4). 

35 


9 


Column 

Code 

Routine/ 

Format 


Explanation 

20 

FLG20 

PARTI /II 

All bodies are surface vorticity bodies 

21 

FLG21* 

II 

Extra crossflow. 

22 

FLG22 

II 

Generated boundary conditions. 

23 

FLG23 

II 

Ring wing option. 

24 

FLG24 

II 

Blank. 

28 

IPUVEL 

II 

Punched output. 

29-30 

NIN 

12 

Unit number for input coordinates 
(defaul t=05) . 

31 

IPRINl 

11 

Matrix print flag. 

32 

IPRIN2 

II 

Matrix-assembly coefficient print flag. 

33 

IPRIN3 

II 

Very detailed matrix construction 
print flag. 

34 

IRAKF 

II 

Automatic rake generation flag (see 
also cards 8 and 9) . 

35 

I SAVE 

11 

Blank. 

1-10 

CHORD 

PARTI/ 

FIO.O 

Reference chord length (default=l .0) . 

11-20 

MN 

II 

Mach number for Goethert correction 
(0.0 implies incompressible). 

21-30 

TCNST 

II 


Blank. 

31-40 

EPSLON 

fl 

J 


1 

IGEOMF 

BASIC!/ 

11 

= 

=0, curved elements 
=1 , flat elements . 

2 

ISIGF 

II 

= 

=0, parabolic a 
=1 , 1 inear a 

-2, constant a (on each element). 

3 

ICURVN 

II 

z 

=0, internally calculated element 
curvatures 

=1 , input curvature (see card 7). 



Card 



Routine/ 



No. 

Column 

Code 

Format 

Explanation 

# 

4 

4 

NONEWF 

BASICl/ 

=0, use the newest formulae 





11 

=1, use the old formulae (implies 






flat elements and constant a). 

% 


5 

I FORM! 

II 

Input format flag (see card 6). 



6-10 

NN 

BASICl/ 

Number of defining end points for 





15 

the body. 



11-20 

MX 

BASIC!/ 

FIO.O 

x-multiplier value (default=l .0) . 

• 


21-30 

MY 

II 

y-multiplier value (default=l .0) . 



31-40 

THETA 

It 

Coordinate rotation value (degrees, 
measured about ~i axis). 

• 


41-50 

ADDX 

II 

x-increment (to be applied to all 
input x-coordi nates for this body). 



51-60 

ADDY 

II 

y-increment (to be applied to all 
input y-coordinates for this body). 

• 

5 

6-10 

BDN 

BASICl/ 

Body number (sequential for bodies, 





15 

zero for off -body points). 



16-20 

SUBKS 

II 

Subcase flag. 



26-30 

NLF 

ff 

Nonlifting flag (for combination 

• 





cases only) . 



31-40 

XE 

BASICl/ 

Semi -major axis for ellipse cases 





FIO.O 

(input if FL607/0) . 



41-50 

YE 

It 

Semi -mi nor axis for ellipse cases 

• 





(input if FLG07/0). 


6 

1-60 

UTXl(I) 

BASIC!/ 

6E13.4 

x-coordi nates . 


6' 

1-60 

UTYl(I) 

It 

y-coordinates. 

i 






6a 

1-10 

UTXl(I) 

BASICl/ 

F10.5 

x-coordi nates . 



11-20 

UTYl(I) 

It 

y-coordinates . 

» 

6b 

1-10 

UTXl(I) 

BASICl/ 

FIO.O 

x-coordinates . 




37 




Card 



Routine/ 



No. 

Column 

Code 

Format 

Explanation 

• 

6b 

21-30 

UTYl(I) 

BASICI/ 

FIO.O 

y-coordi nates . 


7 

1-60 

HCURV(II) 

BASICI/ 

6E13.4 

Curvature values for the NN-1 elements 
which constitute this body. 

• 






8 

6-10 

NN 

BASICI/ 

Number of automatically generated 





15 

mass flow rakes. 


9 

1-10 

XI 

I RAKES/ 
F10.5 

x-coordinate of start of the rake. 

• 


11-20 

Y1 

II 

y-coordinate of start of the rake. 



21-30 

X2 

It 

x-coordinate of end of the rake. 



31-40 

Y2 

II 

y-coordinate of end of the rake. 

• 


41 -45 

N 

IRAKES/ 

Number of intervals to be used in the 





15 

rake (note 4<N<200 and N must be an 
even integer).” 


10 

6-10 

NUN 

BASIC2/ 

Flow identification number. 

• 




15 




16-20 

MSF 

II 

=0, axi symmetric onset flow 
=1, crossflow onset flow 
=2, both 0 and 1 . 



21-30 

TYPE 

II 

= +1.0, velocity will be input in 






x,y component form 






= 0., velocity will be input in 






normal tangential form 






=-1.0, automatic generation of flow 






due to rotation about the 

• 





z-axis (for crossflow only). 



31-40 

FG 

l( 

Flow generator constant. 


11 

1-60 

NG(I) 

BASIC2/ 

6F10.0 

X or normal component velocity. 

• 

ir 

1-60 

TG(I) 

II 

y or tangential component velocity. 




38 


3.2.3 COMBYN Program Input Instructions 


The input instructions for the COMBYN program are taken from 
Reference 1. It is assumed that the user is familiar with the details 
of the COMBYN program. Data file Tape 07, which is one of the output 
data files of the potential flow program EOD, is an input data file for 
COMBYN. Data stored on Tape 07 are the coordinates and the corresponding 
individual velocity solutions of the flow field for on-body and off- 
body points with format (4E13.6). Additionally, Tape 05 generated 
manually by the user is also an input data file for the COMBYN program. 
The card or tape structure of Tape 05 is as follows. 

Card 

Type Description Subroutine 

1 Title card. READS 

2 Number of on-body and off-body points. " 

3] 

> Input initial conditions of the flow. " 

4j 

5 Control flag. " 

6 Circumferential coordinates. " 

7 Location of control station rake. " 

8 Input for defining zero surface distance. " 

English engineering units are used throughout the program. 

Length, in. 

Velocities, ft/sec 
Angles, deg 2 
Pressures, Ib/ft*^ 

Temperature, °R ^ 

Density, slug/ff^ 

Force, lb 

Weight flow, Ib/sec 
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Col umn 

Code 

Routine/ 

Format 

Explanation 

1-72 

TITLE 

READS/ 

18A4 

Title card. 

1-4 

NT(1) 

READS/ I 4 

Number of on-body points for the 
closed-end solution. 

5-8 

NP(1) 

It 

Total number of off-body points. 

9-12 

NT(2) 

11 

Number of on-body points for the 
open-end solution (eliminate the 
last body) . 

13-16 

NP(2) 

II 

Total number of off-body points. 

17-20 

NID 

It 

Number of EOD I.D. cards. 

21-24 

KSKIP 

II 

=0, for first case of COMBYN 
=1, for successive cases using the 
same EOE output. 

25-28 

lOVHUB 

II 

=0, hub vorticity solution from EOD 
is not read. 

=1, hub vorticity solution from EOD 
is read. 

1-8 

VC 

READS/ 

F8.3 

Average axial velocity at the 
control station. Based on live 


flow area, i.e., the flow area minus 
the area associated with the boundary 
layer displacement thickness. If 
WDOT^O, the program will interpret 
this as a code to ignore the input 
VC and will calculate VC from WDOT. 

(To run a case with VC actually 
equal to zero set WD0T=0.0 and 
VC=0.0). If ICTLPT (card 5) is not 
zero, VC is interpreted as the speci- 
fied pressure ratio (PS/Pj) at a 
"control point" rather than a velocity 
Note : All three inputs, WDOT, ICTLPT, 

and VC, must be nonzero when the 
"control point" calculation is desired 

9-16 VINF " Free-stream velocity. 

17-24 ALFAF " Angle of attack, 0.0 for free-stream 

perpendicular to inlet axis. Note 
that ttp = a - 90° . 

For "control point" cases only, 

ALFAF will be calculated when ALFAF 

is input as =999.0. 
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Card 

No. Column 


Code 


Routine/ 

Format 


Explanation 


m 


i 


I 


3 25-32 TTOTAL READS/ 

F8.3 


33-40 ELND 


41-48 YWING 

49-56 UTIP 


57-64 VA 


65-72 PT 


73-80 CUTOFF 


4 1 -8 PSTAT 

9-16 TSTAT 


17-24 WDOT 


Total temperature, if PSTAT and TSTAT 
are read in (to be explained later), 
the program will calculate TTOTAL. 

If TT0TAL=0 and PSTAT and TSTAT=0, 
then TT0TAL=518.67 will be used. 

ELND is the arbitrary length used 
for scaling or normalizing. Refer to 
KND input, card 5. See also CUTOFF 
input below. 

Upper limit of integration for surface 
forces (used in subroutine INFRCE). 

Rotor tip speed. Need not be input 
unless relative rotor inlet quantities 
are desired. (See COMBYN OUTPUT.) 

Bulk velocity at control station, 
i.e., average inlet axial velocity 
based on geometric area. If VA=0.0, 
the program will interpret this as a 
code and set VA=VC. 

Total pressure. If PT=0.0 and PSTAT= 
0.0, the program will set PT=2116. 

If PT=0.0 and PSTAT/0.0, PT is 
calculated. 

If CUTOFF/0, the pressure ratio P^/Pf 
on the shroud will be plotted (on^ 
Calcomp) against dimensionless surface 
distance S/ELND starting at X=XRI and 
proceeding in both directions along 
the surface for a distance of S=CUT0FF 
Length of plot in paper inches is 10 
(CUTOFF/ELND) . There is one plot for 
each circumferential angle THETA. 

Static pressure. 

Static temperature. (If PSTAT and 
TSTAT are not 0.0, total pressure (PT) 
and total temperature (TTOTAL) will be 
calculated using PSTAT and TSTAT. 

Weight flow - required unless VC/O 
and concurrently ICTLPT=0. 
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Routine/ 

Column Code Format 


Explanation 


25-32 DELQ READS/ 

F8.3 


33-40 VPERIN 

41 -48 MC 

49-56 MINF 


Increment in mass flow fraction for 
spacing for calculated streamlines. 
NOTE : Default Value = 0.1, if DELQ 
is input as 0. 


} Blank. 


1 -4 NTHETA 


5-8 NCLO 

9-12 NCHI 

13-16 NX 


READS/14 Number of THETA 's where THETA is the 
circumferential coordinate. If 
NTHETA=0, one THETA (THETA<70°) will 
be read in and used as the“initial 
angle for the start of three- 
dimensional, on-body streamlines. 

For this option, THETA will vary as 
the streamline is followed up the 
shroud instead of remaining a 
constant on one meridian. NOTE : 

No INFRCE (force) calculations or 
pressure plots can be requested when 
NTHETA=0. 

" One rake must be chosen as the 

control station. NCLO is the number 
of the first point on this rake. 

" The number of the last point on the 

control station rake. 

" If NX=-1 , inlet total force calcula- 

tions are obtained (subroutine INFRCE) 
If NX=+1, a supersonic velocity correc 
tion is activated. At those on-body 
points where local supersonic flow is 
detected, velocities and pressure 
ratios are readjusted based on local 
Mach numbers and the rate of change 
of the local velocities. (Since off- 
body points having supersonic velocity 
are not corrected, there will be an 
inconsistency between the corrected 
on-body points and adjacent off -body 
points.) 



Card 

No. 

Column 

Code 

Routine/ 

Format 

Explanation 

• 

5 

17-20 

KND 

READS/14 

Flag for scaling variables before 
velocity and pressure calculations 
and also for nondimensional izing 
after calculations and just before 
printout: 

• 





Scaling: All input lengths and 

coordinates are divided by ELND 
immediately after being input, and 
WDOT is set to WD0T/ELND2. If 
KND = -1, ELND = YTESTS 

= 0, ELND = 1.0 (no scaling) 

= +1, ELND = YTESTS - YTESTH 
= +2, ELND = the read-in value 
from card 3. 

• 





Nondimensionalizing: If KND=8, the 

surface distance, S, will be divided 
by the read-in ELND just prior to 
printout. If KND=9, the on-body X and 
Y coordinates will be divided by the 
read-in ELND just prior to printout. 
NOTE: If CUTOFF is nonzero, surface 

distance will automatically be 
normalized by ELND before printout. 

% 


21-24 

NOTHET 

II 

If = 1 , WDOT and VC will be left 
constant, as input, for all values 
of THETA (neglecting crossflow term). 
If = 0/b1ank, WDOT and VC will be 
corrected for crossflow and will vary 
with THETA. 

# 


25-28 

ICTLPT 

II 

Index number (from EOD output) of 
the desired "control point" where a 
known pressure ratio is to be input 
in lieu of a control station velocity. 
See VC. 

• 


29-32 

ISWIRL 

ii 

(Required when NTHETA=0). Index 
number of point on shroud where 
three-dimensional, on-body stream- 
line calculation will begin; preferably 
near the fan face. 


6 

1-80 

THETA(I) 

READS/ 

10F8.3 

Ci rcumferential coordinate in degrees 
(number of THETA 's read in depends 
on NTHETA, NTHETA<10). 

i 

7 

1-10 

XTEST 

READS/ 

F10.2 

Axial location of control station 
rake. Must be compatible with NCLO 
and NCHI. 
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Column 

Code 

Routine/ 

Format 


Explanation 

11-20 

YTESTH 

READS/ 

F10.2 

Y on the 

hub at XTEST (control station) 

21-30 

YTESTS 

It 

Y on the 
station) . 

shroud at XTEST (control 

1-10 

XRI 

II 

Value of 
distance 

X at which the surface 
is zero. 

11-20 

YRIHUB 

II 

Y on the 

hub at XRI . 

21-30 

YRISHR 

II 

Y on the 

shroud at XRI. 

31-34 

NHUBMX 

READS/14 

The number of the last point on the 


hub (this can be found in the printed 
output of SCIRCL) . 


3.2.4 Particle Trajectory Program Input Instructions 


m 


# 


There are five input data files used in the trajectory program. 

They are described in the following: 

1. Tape 05 stores the airfoil geometry data. This is the 
same file as the input file Tape 05 of the potential 
flow program, Section 3.2.2. 

2. Tape 21 stores the data for calculating the flow 
velocities. This is the same file as the output file 
Tape 21 of the potential flow program. 

3. Tape 15, which is one of the output data files of 
COMBYN program, stores the coefficients for calcu- 
lating the combination solution. 

4. Tape 02 stores the initial input data for the trajectory 
program. Details are described in the following. The 
data must be in MKS units with the exception of 
particle diameter, which is input as microns. 

The input deck structure for the particle trajectory program is as 

follows. 

Card 

Im. 

1 

2 

3 

4 

5 ' 

6 

• 

7 

8 

J 

9 

10 ' 

' ■ 

10 

/ 


Description 

Number of bodies; degree of governing equations 
Flow field control flags 
Trajectory control flags 
Initial conditions of flow field 

Initial conditions of particles 

Number of size increments 

This card is input only when NSI>1; number of 10 cards = 
NSI 
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Card 

Type Description 


11 

12 

13 

14 

15 1 

15 

16 ■ 

17 

18 
19 . 


Input for polynomial least square fit 
TITLE of experimental data 
Number of experimental data 

Number of 13 cards = NASA 


Input for plotting 


For more than one body, cards 11 through 19 are cycled NBDY times. 


Card Routine/ 

No. Column Code Format Explanation 


1 6-10 NBDY 


16-20 NPL 

26-30 NSEAR 


36-40 NEQ 


Main/15 Number of bodies. 

=1, One-body case (airfoil without 

walls) 

=2, Two-body case (inlet without hub) 
=3, Three-body case 

i) Airfoil with walls (LNTL=1) 
ii) Inlet with hub (LNTL/1) 

Axi symmetric inlet case, NBDY=3. 

" Number of particle trajectories to 

be computed. 

" Maximum number of loops allowed in 

the search for the upper and lower 
impingement limits for the case 
LIM=1 (card 3). 

" Number of equations to be solved. 

NEQ=6 for a lifting, rotating 
particle; NEQ=4 for a spherical 
particle undergoing drag only 
(Section 2.0). 


Card Routine/ 

No. Column Code Format 


Explanation 


2 6-10 LEQM 


16-20 LSYM 


26-30 LRANG 


3 6-10 LIM 


16-20 LOPT 


26-30 LPLOT 


36-40 LTNL 


46-50 LXOR 


Main/15 =1, The initial particle velocity 
is equal to the flow at the 
initial particle location. 

, The initial velocity conditions 
input on cards 4 and 6 are 
= used as the initial particle 
velocity. 

" Symmetric flow field flag 

=0, Unsymmetric flow field (general 
case). 

=1, Symmetric flow field (only 
half plane is computed). 

Axi symmetric case LSYM=0. 

" /O, Locates approximate values of 

yQM and yoj;, (see Section 
3.1.6). 

" /O, Calculates surface impingement 

limits you and yg£. 

=0, Calculates single particle 
trajectories; LRANG=0 (see 
Section 3.1.6). 

" /O, Stores details of particle 

trajectories on data Tape 04. 

" /O, Executes subroutine for plotting 

local collection efficiency, 6, 
versus surface distance, s/t^. 

" =1, Computes flow over airfoil with 

walls (i.e., wind tunnel simulation) 
/I, Computes flow over airfoil without 
walls and inlet case. 

Axisymmetric case, LTNL/1 . 

" Control flag for adjusting the 

upstream x-coordi nates, XORC, for 
upstream position of particle 
release (see Section 3.1.3). 

=0, The upstream x-coordi nates of 
particle release is not adjusted. 
Particles are released from XORC 
input on card 6 by user. 

=1, Particles are released from the 
x-coordi nate position adjusted 
by the criteria: 

|1 - (W/W )| < EPSX. 
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Column 

Code 

Routi ne/ 
Format 

Explanation 

6-15 

G 

Main/ 

FIO.O 

Acceleration of gravity, 
m/$2 . 

21-30 

YINF 

II 

Velocity of free-stream, 
m/s. 

6-15 

DP 

II 

Particle volume median 
diameter in microns (10"^ m) 

21-30 

RL 

II 

Reference length, e.g., 
chord length of the airfoil, 

m. 

36-45 

TIMSTP 

II 

Initial value of the time 
step used in the Adams- 
Moulton predictor-corrector 
method (Section 2.3). 

51-60 

YLIM 

II 

Accuracy criteria for 
computing the surface 
impingement limits (Section 
3.1.6) YLIM=10-5 m is 
recommended. 

66-75 

CF 

II 

Cunningham correction 
factor. Use if volume 
median diameter DP<10y; 
otherwise, CF=1 . 

6-15 

ATK 

II 

Initial value of the 
angle, a. Figure 2.1 (deg). 

21-30 

PIT 

II 

Initial value of the angle, 
0, Figure 2.1 (deg). 

36-45 

PITDOT 

II 

Initial value of 0 (deg/sec) 

51-60 

XORC 

II 

Upstream x-coordinate posi- 
tion of particle release, 
x/i^ (LX0R=0). 

66-75 

YORC 

II 

Upstream y-coordinate 
position of particle 
release, y/ic (NPL=1 , 
LRANG=0, and LIM=0). 

6-15 

XREAR 

II 

Maximum downstream value 
of x/Iq for which the 
flow field is computed. 


Card 

No. 

Col umn 

Code 

Routine/ 

Format 

Explanation 

7 

21-30 

YYLO 

Main/ 

FIO.O 

Minimum value of y/J.^ for 
which the flow field is 
computed . 


36-45 

YYUP 

II 

Maximum value of for 

which the flow field is 
computed . 


51-60 

YMAX 

II 

Initial guess of the upper 
limit you/J^c ) 

(Section 3.1.6). 


66-75 

YMIN 

II 

Initial guess of the lower 
1 imit (LIM=1 ) 

(Section 3.1.6). 

8 

6-15 

ADEPS 

II 

Convergence criteria of 
the Adams-Moul ton predictor- 
corrector method. 

ADEPS=0.001 is recommended. 


21-30 

ANLFLO 

II 

Angle of flow direction 
relative to airfoil chord 
line measured from positive 
direction of x-axis. 


51-60 

EPSX 

II 

Accuracy criteria for the 
case LX0R=1 (see card 3); 
EPSX=0.001 is recommended. 

9 

6-10 

NSI 

Main/15 

Number of droplet size 
i ncrements . 

10 

This card 

is input 

only when NSI>1 . 



1-10 

PLWC 

Main/ 

FIO.O 

Percentage of liquid water 
content for each drop size. 


11-20 

DPD 

II 

Distribution of the particle 
diameter ratio to volume 
median diameter. 


21-30 

CFP 

II 

Cunningham correction factor 
corresponding to DPD (see 
Table 6.2). If DP*0PD>10u, 
CFP=1. 

Number 

of 10 cards 

= NSI. 
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Card 

No. 

Column 

Code 

• 

11 

1-80 

TITLE 


12 

6-10 

NCOEF 



21-25 

NPTS 



36-40 

NS 


13 

1-80 

TITLE 


14 

1-5 

NASA 


15 

1-10 

SPOIL 



11-20 

EFC 


16 

1-80 

XTITL 


17 

1-80 

YTITL 


18 

6-15 

YHEIT 



21-30 

YMAX 



36-45 

XMIN 


Routi ne/ 
Format 

Explanation 

Mai n/ 
20A4 

Initial input data for 
local collection efficiency 
calculations. 

Main/15 

Number of coefficients of 
polynomial curve fit for 
calculating B. The order 
of the polynomial is NCOEF-1 
(Section 3.1.7). 

II 

Number of data points used 
in curve fitting for com- 
puting local collection 
efficiency (Section 3.1.7). 

II 

=0, Segment curve fit data 
=1, Total curve fit (see 
Section 3.1.7). 

Main/ 

20A4 

Experimental values of local 
collection efficiency. 

Main/ 1 5 

Number of experimental data 
points input for comparison 
with computed results. 

Main/ 

F10.5 

Surface distance, s/Iq, 
along the body at the 
location of the measured 
value of local collection 
efficiency. 

II 

Experimental values of 
local collection efficiency 
at positions SPOIL. 

APLOT/ 

20A4 

Title of the x-axis. 

II 

Title of the y-axis. 

APLOT/ 

FIO.O 

Length of y-axis for graph 
to be plotted (inches). 

II 

Highest value of datum 
point on y-axis (inches). 

II 

Far right value of datum 
point on x-axis (inches). 
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Card Routine/ 


No. 

Column 

Code 

Format 

Explanation 

19 

6-15 

XLENG 

APLOT/ 

FIO.O 

Length of x-axis for graph 
to be plotted (inches). 


21-30 

XMAX 

tl 

Far left value of datum 
point on x-axis (inches). 


36-45 

XMIN 

APLOT/ 

FIO.O 

Far right value of datum 
point on x-axis (inches). 


Cards 11 through 19 are cycled NBDY times. 

3.3 Description of Program Output 

The output files of the programs are divided into three categories: 
(1) output of the geometry generation program, (2) output of the axisym- 
metric potential flow program, and (3) output of the particle trajectory 
program. For the first two programs only the output files related to 
the present particle trajectory program are described. 

3.3.1 Geometry Generation Program Output Description 

The geometry generation program is named SCIRCL. The output data 
file Tape 17 of SCIRCL is used as the input data file Tape 05 of the 
axi symmetric potential flow program and of the particle trajectory 
program. Tape 17 stores the airfoil geometry data. In general, Tape 17 
is the same for both the EOD and the particle trajectory program; 
however, in some cases it is modified when used in the trajectory 
program as described in Section 3.2.2 and Section 5.3. The format of 
Tape 05 is the same as the input data described in Section 3.2.2. 

3.3.2 Axisymmetric Potential Flow Program Output Description 

The output data file Tape 21 of the potential flow program is used 
as the input data file Tape 21 for the particle trajectory program. 

Tape 21 records all the necessary information including the source, 
sink, and/or vorticity distributions along the body surfaces for calcu- 
lating the flow field about the bodies. Note that Tape 21 will not be 
generated by the axisymmetric potential flow program unless at least one 
value has been assigned to NRAKE during input to the SCIRCL program. 


Thus, file Tape 05 must contain one or more off-body points (Section 
3.2.2). 

The output data file Tape 07 of the EOD program is used as the 
input data file Tape 07 for the COMBYN program. Tape 07 records the 
coordinates and the corresponding individual solutions for on-body and 
off- body points. 

3.3.3 COMBYN Program Output 

The output data file Tape 15 of the COMBYN program is used as the 
input data file Tape 15 for the particle trajectory program. Data 
stored on Tape 15 are the coefficients for evaluating the combination 
solution. 

3.3.4 Trajectory Program Output 

A data file Tape 04 is generated if LOPTj^O and contains the values 
of X, y, X, y, 6, and e for each time step of the particle trajectories. 
Values of the flow field velocities at the particle locations along the 
particle trajectory are also contained on this file. 

Data file Tape 08 is also generated if LOPTj^O and contains the 
particle locations along the computed trajectories for plotting purposes. 
The data stored on Tape 08 and Tape 04 are nondimensional . Because Tape 
04 and Tape 08 are used to record the information at each time step, the 
program will need much more storage than that generally needed with the 
option L0PT=0. LOPT/0 should be used only if a few particle trajec- 
tories are to be calculated. 

Additionally, information on particle initial position, impingement 
limits, surface impingement point coordinates, and surface distance are 
recorded on data file Tape 06 for printout. The computed local collec- 
tion efficiency, e, is also recorded on Tape 06. The value of B versus 
surface distance is recorded on Tape 09 for plotting purposes. 

Tape 09 is unformatted. The format for Tape 06 will be explained 
in Section 5.0 where test cases are documented. The format for Tape 04 
and Tape 08 is illustrated as follows. 


Format for Tape 04 (F0RMAT(1QE10.3)) 


Column 

Code 

Explanation 

1-10 

X 

Time 

11-20 

El 

x-coordinate of the particle, x/ji^ 

21-30 

E2 

y-coordinate of the particle, 

31-40 

YLAST(5) 

Pitch angle of the particle in radians 

41-50 

E3 

x-component velocity of the particle 

51-60 

E4 

y-component velocity of the particle 

61-70 

YLAST(6) 

z-component angular velocity of the particle 

71-80 

E5 

Velocity of the particle relative to the flow 
field 

81-90 

W1 

x-component velocity of the flow field 

91-100 

W2 

y-component velocity of the flow field 

Format 

for Tape 08 (F0RMAT(3F10.5) ) 

Column 

Code 

Explanation 

1-10 

XP 

x-coordinate of the particle, x/t^ 

11-20 

YPL 

y-coordinate of the particle, x/t^ 

21-30 

SW(IW) 

The parameter is used to indicate the end of 
each trajectory, i.e., if SW?^88.8888 the 
trajectory is terminated. 


4.0 PARTICLE TRAJECTORY COMPUTER PROGRAM CAPABILITIES 
AND FUNCTION DESCRIPTIONS 


Details of the particle trajectory computer program are presented 
in this section. Each of the subroutines is described individually in 
Section 4.1. A list of FORTRAN variable names is given for each sub- 
routine, and flowcharts are provided to illustrate the order of calcula- 
tion. Error messages and typical corrective measures are given in 
Section 4.2. 



VINE Free-stream velocity, W 

YP y/£^ coordinate of current particle location 

YHIT Initial coordinate of a particle which impacts body 

YMIS Initial y/Jl coordinate of a particle which moves away from 

the body ^ 

SL Lower surface impingement limit, 
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su 


Upper surface impingement limit, 

ZL Lower limit, corresponding to SL 

ZU Upper limit, corresponding to SU 

See Figure 4.1 


Flowchart: 



Read Cards 1 
10 of TAPE 02 



Search for 
Lower 

Impingement 
Limi t 


Search for 
^ Upper 
Impingement 
Limi t 

as 


YES 


IW=0 

IC=0 

IT0T=NPL 





Read Cards 1 1 - 
1 5 of Tape 02 



Read 

Cards 

16 

- 19 


( Call 

APLOT 


Call 

TRAJCT 



• Set 

Equil ibrium 
Initial 
Conditions 


Figure 4.1 Flowchart of main program. 
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4.1.2 Subroutine ADAMS 


Objective: 

Use Adams -Moulton, variable time 
method to solve the equations of 

step, predictor-corrector 
particle motion. 

Variables : 



ypred 

Adams-Moul ton predictor 


YCORR 

Adams-Moul ton corrector 


ylast 

Value of function y (see Section 

4.1) at previous time step 

TSP 

Initial time step increment 


X 

Time 


Flowchart: 

See Figure 4.2 





ADAMS 


V 



Figure 4.2 Flowchart of subroutine ADAMS. 
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4.1.3 Subroutine RANGE 


Objective; 

Flowchart: 


Locate an approximate range of which contains the 
limits and 

See Figure 4.3 


Y0(1 j=YMX 
Y0(2)=YMN 
DY=Y0{1)-Y0{2) 
IAC=0 
IW=0 
IPR0T=0 



^ / 

- 

IW=IW+1 

IPR0T=IPR0T+1 







Set Initial 
.Conditions of 
the Particle 
for LEQM=1 Case 


WM 



r 

YL(IW)=YP-YRER 


_ 



Yes 



Y0(2)=Y0(2)-DY 

IAC=0 

IW=1 



YMX=Y0(1) 

YMN = 

Y0(2) 1 

1 

r 


Figure 4.3 Flowchart of subroutine RANRE. 
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c 


RETURN 


) 














4.1.4 Subroutine F 


• 

Objective: 

Supply functional form of equations 
motion. 

governing the particle 


Variables : 



• 

T 

X 

Time 

Values of the functions (correspond 

to Y in Section 4.1 ) 


XDOT 

Derivatives of x with respect to time (i.e., x, y, e, x, y, 
6, see Section 2.2) 

• 

CL 

Coefficient of lift force 



CD 

Coefficient of drag force 



CM 

Coefficient of moment 


• 

4.1.5 Subroutine COEFF 



Objective: 

Supply aerodynamic coefficients CL, 

CD, and CM. 

• 

4.1.6 Subroutine MODE 



Objective: Determine whether the particle impacts the body. Also, 

terminate the trajectory calculation for particles which 
move away from the bodies. 

Flowchart: See Figure 4.4 
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Figure 4.4 Flowchart of subroutine .MODE 











4.1.7 Subroutine READIN 


Objective; 
FI owchart; 


Read Tape 21 and Tape 05. 
See Figure 4.5 



Figure 4.5 Flowchart of subroutine READIN. 







4.1.8 Subroutine VELCTY 


Objective: 


Variabl es : 

X 

Y 

VXC 

VYC 


Compute flow velocities for a given position (x,y). 
Subroutine VELCTY calls subroutines MATRIX, AXIS, CROSS, 
SINTP, and VBARIT directly for each velocity calculation 
Subroutines SORTXY, ARSIN, ELINT3, ELIP, ELLC, HLAMB, 
INEL, PARAB, PARAB2, QC, and SIMSON are called from EOD 
and COMBYN only once. 


Dimensionless x coordinate, x/s,^, of particle position 
Dimensionless y coordinate, y/«.j.> of particle position 
Dimensionless x component velocity, W /W , of airflow 

X 

Dimensionless y component velocity, of airflow 


Flowchart: See Figure 4.6 







Figure 4.6 Flowchart of subroutine VELCTY. 
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4.1.9 Subroutine EFFICY 


Objective: 

Variables : 

NCOEF 

NPTS 

NS 

YP 

XIN 


Calculate the local collection efficiency, b, using both a 
linear approximation and a polynomial curve fit interpola- 
tion technique. 


Number of coefficients for polynomial curve fit 

Number of data points used for each curve fitting procedure 

Do-loop index (see flowchart. Figure 4.7) 

Collection efficiency, B = -dy^/ds 
Value of sCy^) 

See Figure 4.7 


FI owchart: 


Call LINEAr\ 
(linear ) 
interpolation)/ 


Call TERP 
(polynomial 
fitting) 


I 


Record Data on 
Tape 09 
(for plotting) 



Figure 4.7 


Flowchart of subroutine EFFICY. 
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4.1.10 Subroutine EF 

Objective; Calculate the local collection efficiency, 6, for multi-size 
particle distribution NSI>1 case using a polynomial curve 
fit interpolation technique. 

4.1.11 Subroutine LINEAR 

Objective; Compute the local collection efficiency by a linear 
approximation. 

Variables ; 

X Value of s 

Y Value of y^ 

YP Local collection efficiency, -dy^/ds 

4.1.12 Subroutine TERP 

Objective; Curve fit the {s,y^} data to a polynomial function. 

Variables ; 

XIN x(s) coordinate at position (x,y) 

YOUT y(s) coordinate at (x,y) obtained by the polynomial data fit 

YPOUT Derivative of y with respect to x, dy/dx (dy^/ds) 

XA Data set of x-coordinate 

YA Data set of y-coordinate 

N1 Index of the first data point used for data fit 

N2 Index of the last data point used for data fit 

NCOEF Number of coefficients in the polynomial function 

NPTS Number of data points used in each curve fit procedure 

4.1.13 Subroutine CHOLES 

Objective; Matrix solver called by subroutine TERP. 
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4.1.14 Subroutines APLOT and TRAJCT 


Objective: Plot routines for 6 and for particle trajectories to be 

furnished by intended users. 

4. 2 Error Messages 

Table 4.1 lists the error messages and recommended corrective 
actions to be taken for the computer program. The particular subrouti 
in which the error message occurs is also listed in the table. The 
table is essentially self descriptive and needs no further discussion. 


TABLE 4.1 Error Messages and Corrective Actions. 


Error Message Generated by 
the Program 

Problem and Corrective Action 

Error 

Location 

The number of y- 

Corrective Action: Check 

the 

Subroutine 

coordinates read is not 

input data file. Tape 05, 

for 

READIN 

equal to the number of 

equal number of x- and y- 



x-coordinates read. 

coordinates . 




Time-halving loops exceed Problem : The Adams-Moul ton Subroutine 

100. method of solution is not ADAMS 

converging to the desired 
accuracy at the current time 
step. The problem itself 
sets NL00P=100. The value is 
large enough for all cases 
investigated in this report 
where TIMSTP is set equal to 
10 - 4 . 

Corrective Action : 

1) Reduce the magnitude of 
TIMSTP. 

2) Increase the size of 
ADEPS. 

3) Check for severe gradients 
or error in the flow field. 


IP0RT.GE.31 in SUB.RANGE-- 
program stopped. 


Problem : Initially guessed 

values of YMAX and YMIN are 
too far from the body or 
AY = YMAX - YMIN is too 
smal 1 . 

Corrective Action : 

1 ) If particles pass below 
the body, increase the 
input values of YMAX and 
YMIN. 

2) If particles pass above 
the body, decrease the 
input values of YMAX and 
YMIN. 

3) Increase magnitude of aY. 


Subroutine 

RANGE 



TABLE 4.1 (continued). 


Error Message Generated by 
the Program 

Problem and Corrective Action 

Error 

Location 

FORMAT 172: Number of 

Problem: Number of particle 

MAIN 

particle trajectories 

trajectories which hit body 


exceeds NSEAR, program 

exceeds the input value of 


stopped. 

NSEAR (Card 1 on Tape 02). 

Corrective Action: 

1) Increase the input value 
of NSEAR. NSEAR=50 is 
recommended. 

2) Increase the input value 
of YLIM (Card 5 on Tape 
02). 

3) Decrease the input value 
of NPL (Card 1 on Tape 
02). 



FORMAT 747: Number of 

trajectories is more than 
100, program stopped. 


Problem : Number of trajec- MAIN 

tories calculated exceeds 
dimension state. Either the 
program is having difficulty 
finding limits or specified 
value of NPL is too large. 

Corrective Action : 

1) Refer to the actions in 
error message #4. 

2) Adjust input values of 
YMAX and YMIN (Card 7 
on Tape 02--refer to 
error message #3) . 


Particle released between 
the region ZL-ZU crosses 
upper or lower impingement 
limit trajectory, program 
stopped. 


Corrective Action : Check 
the flow field around the 
body and increase absolute 
value of pseudo-surface DX 
(Card 2 on Tape 05). 


MAIN 
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5.0 TEST CASES FOR PARTICLE TRAJECTORY CALCULATIONS 


Two test cases for water droplet trajectory calculations are given 
in this section. Excessive CPU time is required with the EOD program, 
even with the significant improvements to accelerate the computation as 
described in Section 6.0. For this reason calculation of local collec- 
tion coefficients is impractical. The calculation of sufficient trajec- 
tories, however, has been carried out for the test cases to fully verify 
the trajectory program. The test cases illustrated are: 

1. Axisymmetric inlet at 0° angle of attack with M = 0.4. 

2. Axisymmetric inlet at 30° angle of attack with M = 0.4. 

An interpolation scheme was tested but found to have no advantage 
in terms of computational time for the specified test cases and, therefore, 
is not recommended (see Section 6.0 for details). 

The input/output step-by-step procedures for applying the program 
to each test case is described. Details of test cases 1 and 2 are 
described in Sections 5.1 and 5.2, respectively. To fully understand 
the input/output printouts presented in these sections, the user should 
coordinate the results with the card structure and input instructions 
given in Section 3.2. 

A source listing of the particle trajectory computer program is 
given in the appendix. 

5. 1 Axisymmetric Inlet at 0° Angle of Attack with Mach Number = 0.4 

The geometry of the axisymmetric inlet is given in Figure 5.1. The 

input data for generating the geometry using the program SCIRCL is 

listed in Table 5.1. The output data file Tape 17 of SCIRCL (Table 5.2) 

is used directly as the input data file for the axisymmetric potential 

flow program (Section 3.2.2). The output data file Tape 07 of the 

potential flow program and Tape 05 (Table 5.3) created manually by the 

(Section 3.2.3) are used as the input data files for the COMBYN 
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user 


Figure 5.1 Geometry of axisymmetric inlet normalized by inlet diameter 
(1 1 .68 units) . 


TABLE 5.1 Test Cases Computed with Particle Trajectory Program, 


:£i£'T CASE RELEASE 2.0 


r G L 1 

EDO 


2.00 

0.02 0.08 

1 

* .0425 

0.0 0 « Si i i 5 

10 

j . 0 

0.0 0.4709 

10 

j . 0 4 25 

0.0 0.0 

10 

0 . T. **• - 

0 .0 0.0 

io 

0 . i .? c 4 

0 .0 0.0 

10 

J . i 4 5 4 

0.0 0.0 

10 

0 . 

0.0 0.0 

10 

j . 3 1 tr .S 

0 . 0 0 . 0 

\0 

5 . 997.2 

0.0 0.0 

10 


2.0 2.0 


■ ■ * cSZc? 

0 .685 87 


" 1 . ') 0 

0 . 00 


1 . 0 0 
j t 9 9 A09 

3 . 0 S 2 1 9 


0 . 205-.O 

0 . 2C54S 


. 0 

6 . 0 


TToacis 

0 . 99409 


• ; , 5 1 7 0 

0.51370 


'• .: » \j 0 

0.99409 

0 . 1 4538 

0.51370 

0.51370 

0 . 42731 

; 000 . 0 

2 . 2 . 0 


' . » j i c j 

0 . 1 4538 


2' . *12 7 21 

0.42 /' J 1 


1000 . 
. 0 0 

0 . 00 

0 . 02209 


0 . 500 

0 . 52389 

■ . 

» *.27098 

0.99409 


• ^ 2 -1 h 

0 .55248 


1 . •} 0 
. , , c- 4 0 9 

3 .08219 


0 . j t) ^ 4 S 

0.55248 



0.00 


0 .06627 
0 .537A1 


0.99409 
0 .20548 


0 .04281 
0.42731 

0.00 

O.S 

0 . 22098 
0.55248 


3 .08219 
0 .20548 


0.00 

0.6 

0.34247 

0.55248 
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TABLE 5.2 Input Data File for Geometry Generation Program SCIRCL 
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0.51270E+00 0.51370E+00 0.51370E+00 
0.51370E+00 0.51370E+00 0.51370E+00 
0.51370E+00 0.51370EE00 0.31370E+00 
0 .51370E + 00 O.5137OEt00 0.513 70E + 00 
0.51370E+00 0.51370E+00 0.51370E+00 
0.51288E+00 0.51227E+00 0.51150E+00 
0 . i0698E+00 0 . 5O533ErO0 0.5039EE+00 
0.49663E't00 0.49457E+00 0.49243E+00 
0.48313E+00 0.48073E+00 0.47R24E+00 
0.46798E+00 0.46541E+00 0.46287E+00 
0.45313E+00 0.45084E+00 0.44861E+00 
0.44033E+00 0.43874EE00 0.43705E+00 
0.43150E+00 0.43043E+00 0.42951E+00 
0.42740E+00 0.42731E+00 0.42783E+00 
0.43354E+00 0.44307E+00 0.44832E+00 
0.47523E + 00 O.43323iEfO0 0.49173E + 00 
0.51833E+00 0.52170E+00 0.5248 7 E+00 
0.S3601E+00 0.53350Et00 0.54089E+00 
0 . 54902E+00 0 . 55042E+00 0.5514SE+00 
0.55248E+00 0.55246E+00 0.55248E+00 
0,55243Et00 0.55248EEO0 0.55248E+00 
0.S5248E+00 0.55248Et00 0.5524BEt00 
O.S5248EtOO 0.55248E+00 0.5524HE+00 
0.55248E+00 0.55246Et00 0.55246E+00 
0.55243EfOO 0.55248E+00 0.55248E+00 
0.S524aET00 0.5524SE+00 0.55248E+00 
0.55248E+00 0.55248E+00 0.55243E+00 
O.S5 24£tiTOO 0.55 2 48E + 00 0.5524BE + 00 
0.55248EEOO 0.55248E+00 0.55248E+00 
0.55248E+00 0.55248E+00 0.5524HE+00 
0.55248E+00 0.55248E+00 0.55248E+00 
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TABLE 5.2 (continued) 


0000 90 




EFOCl 


0 

0 

0 


8 F 0 C 1 


-0 . 42800E-01 

-0. 42800E-01 

-0 . 42300E-01 

-0. 42300E-01 

-0.42800E-C’l - 

0 . -i C* 0 0 E •• 1 

-0 . 42300E-01 

-0.42300E-01 

-0.42300E-01 

-0. 42300E-01 

0 , OOOOOEtOO 

O.OOCOOEiOu 

0 . OOOOOE+00 

0 . OOOOOE+00 

0 .OOOOOE+00 

0 . COOOOEtOO 

0 . OOOOOEtOO 

0 . C C> ('■ C C E f 0 •!’ 

0 .OOOOOE+OC 

0 .OOOOOE + 00 

0 . 42300E-01 

0 .42300E-01 

0 .42S0OF.-01 

6 . 43 .” C- C- E - 0 1 

0 . 42800E-01 

0 . 42800E-01 

0 . 42SO0E-O1 

'' . 42800E-01 

0 . 42800E-0 1 

0 • 43 • ' 0 E “ 1 

0.85A00E-01 

0.S3400E-01 

0 .3E.600E-01 

0.85A00E-01 

0.95600E-01 

0 .=5r00£-0> 

0.85&00E-01 

0.85600E--01 

0 .S^.600E-0i 

O.S^600E-01 

0. 12840Ft00 

r* . < " - u - - 1 - . ' 

0 . 1 2840E+00 

0 . 12340E + 00 

0. 12840E + 00 

0 . 12340Et00 

0 . l2840Et00 

0 . 1 Z c 4 0 £ r C 0 

0 . 1 2840E+00 

0 . 1 2340E+00 

0 . 1 4540E + 00 

0 . 14540EtOO 

0 . 1 4340 Ft 00 

Q , 1 4 ' A } £ + J 

0 . 1 4".40E + 00 

0. 14340E+00 

0 . 14540E toe 

0. 14540Et00 

0 . 1 4‘340Ft6c 

C' . i 4 5 4 : E * -3 

0 . E'.UOOEtOO 

0 . 23 lOOE+00 

e .::3iooEtjo 

■J . 2.3iOOEtOO 

0 . 2 3 1 ‘1 0 " t 0 0 

•J . Z 4 i i. \ C 

0.23100E+00 

0 . 231 OOE + 00 

0 . 23100Et00 

■7 . 231OOEtO0 

0 . 31^FOEtOO 

0 . 3 i 6 r ‘J £ : 

0.!16^.0E + 00 

0 . 31640E + 00 

'■ . 31 &60E too 

0 . ilAAOEtOO 

0.31 6F OFt'"-0 

r. ^ T 1 4 4 . ^ r.' l- r. 

0 . 21660E+00 

0. 3l640E + 00 

0.'>'V320r:t'.iO 

O.9?32OEt0O 

0 . r9320£t00 

► “ ’•7 J. w £ r Z 0 

0 . 95320Ef00 

0 .99320F+00 

0 . 99 3 20E too 

0 .99320Et00 

0 , -7 9 1 ~ 0 1 -' 4 0 

r. c 3. r. r. 

0 . OCOOOE+00 

0 . 41833E-01 

0 . 12367EtOO 

0 . 18550Et00 

0 .24733E400 

0 . 3 0 7 i .7 £ * 0 

0 . 37100E+00 

0 . 43283E + 00 

0 . 4?467Et00 

0.55650Et00 

O.OOOOOEtOC 

0 . 5 3 3 3 Z £ - Z 1 

0 . 10464E+00 

0 . 15&97E + 00 

0 .20929EtOO 

0.26l61Et00 

0 .313?3£t00 

Z • 3 ^ c £ T C' . 

0.41S53E+00 

0 . 47090E + 00 

O.OOOOOEtOC 

0.43489E-01 

0 .94979E-01 

■' . 1 *» 'j ^ T 0 C' 

0 . 19396E + 00 

0 . 24244E+00 

0 . 2 9093Et00 

0 . 33942Et00 

0.3S791Et00 

0 . 4364e3£ + C0 

0 . OOOOOE + 00 

0.44534E-01 

0 .93067E-0’ 


^ , 1 9 6 1 3 F 1 0 

r , •> ;7 2 ~ cr • 

0 . 27920E + 00 

0 . 32373E tOC- 

0 .27227E400 

0 . 4 LSSOEtdo 

0 .OOOOOEtOO 

•3 . 4 3 J - E - Z 1 

0. 91217E-01 

0 . 136S3E too 

0 . lo243EtOO 

0 . 2280 4E too 

0 . 2 7 3 6 j F r 0 C 

0 * v’ 1 ' Z r £ 7 '*• Z 

0 ■ 3 6 4 8 7 E 4 0 0 

j . 4 1 0 4 8 E T 0 0 

0 . 0 0 0 0 0 E 1 0 0 

0 . 45.36&E--01 

0 . 9 1 3 3 1 E - 0 1 

0 » 1 5 "T" _ £ r 0 Z 

0 . 1 82i4E + 00 

0. 22833Eft-0 

0.27399EtOO 

0. 31943E+00 

0 . 3 A f: 3 3 E t C 0 

0 . 1 Z’ ‘ * £ T 0 C 

0 . oooocE + ;o 

0 . 45962E-01 

0 .91924E-01 

0 . 137S9Et00 

0 . l3335EtOO' 

Z"' 7 r .• 

0 . 27577E + 03 

0.32173E+00 

0 . 34770EtOO 

0.41346EtO0 

0 - OOOOOEtOO 

r , , 4 .i. 4 ; T ■« - r ^ 

0 . 9 3337E-01 

0 . 1 4008E + 00 

0 . 18477EtOO 

0. 2334 7E too 

0 . 2 £ 0 1 6 £ 1 0 0 

. 3zic?iiC0 

0 - 373^nE + 00 

0 . 42024E + 00 

0 ,22195EtOO 

0.25254EtOO 

0 . 293 17EtOO 

r T t r - c c X r •• 

■' . 344 .)OE too 

2 . 3730 IE too 

0.40562EtOO 

0 . 43623Et00 

0 . 4 A 1 8 3 E 4 0 0 

. 4 - 4 6 E, 1 0 6 


TABLE 5.3 Input Data File Tape 05 for the COMBYN Program, 



TEST CASE RELEASE 2.0 

242 90 

242 90 1 0 1 

0 

434.00 

145.30 -90.000 458. 

180 0.00 

.000 

.000 0.0000 

100 

2 81 

90-1 0 1 0 

0 

.000 

180.000 


0 . 9932 

0.2035 0.5137 


.00 

.00 0.50 

55 


0.000 


. 000 


? 1 4 . j i 0 


computer program. The output data file Tape 17 (renamed as Tape 05) of 
program SCIRCL, the output data file Tape 15 of program COMBYN, and the 
output data file Tape 21 generated by the potential flow program are 
used as the input data files for the trajectory program (Section 3.2.3) . 

In addition to the data files described above, the user must 
manually create data file Tape 02 (Table 5.4) for the trajectory program. 

In data file Tape 02 (Table 5.4), the input value of NPL is 1 and 
the particle initial position is assigned at x = -0.5 and y = -0.7. In 
this case, we only consider one particle trajectory which is released 
from the prescribed position as data file Tape 02 (Card 6). The particle 
hits body 1 at x = 0.01325, y = -0.51912 with surface distance s = 

0.02432 measured from the leading edge point. The particle released 
from X = -0.5 and y = -0.71 will miss the body from the lower side. 

Four particle trajectories are shown in Figure 5.2. 
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TABLE 5.4 


Input Data File Tape 02 for the Particle Trajectory Program 
(0° angle of attack). 


t N 8 D Y 

3 

NFLtt 

1 

NSEAR 

60 

NEQtt 

4 


tLEQM 

1 

LSYM4 

0 

LRANG 

0 




tLIMt 

0 

LOFTt 

0 

LFLOT 

0 

LTNLt 

0 

L 0 r*. i 

to * 

0.0 

UINFt 

44.35 






tUF t 

16.7 

RL t 

1 .00 

TIMSF 

0.0010 

YLlMt 

0.00000 

ICr * 

FRATK 

0.0 

FIT t 

0.0 

FITDT 

0.0 

XOFCt 

-C .5 

Y 0 R C T 

XREAR 

0 . 40 

YLO t 

-0.60 

(UF t 

-0 . 4 

Y MAX t 

r. , -1 

< rt I m 

ADEFS 

0.0010 

ANLFL 

0 . 0 

ANLCR 

0 . 0 

EFSXt 

o"ooi 


NSItt 

1 








INITIAL 

INFUT 

DATA FOR 

LOCAL 

catch 

EFt^ICIENC 

Y CALC' 

L'LAT IONS 



NCOEF A NFTSt 10 1 Hot 0 

EXFERIhENTAL tiATA OF LOCAL CATCH EFFICIEMCY : 


1 

-0.08 0.045 

INITIAL INPUT CiATA FOR LOCAL CATCH EFFICIENCY CALCULATIONS 
NCOEF 4 NFTS» 10 i NS* 0 

EXP-ER IMFNTAL DATA OF LUCAL CATCH EFFICIENCY I 

I 


-0.03 C . 0 4 5 

INITIAL IN'^-Ut' DATA FOR LOCAL CATCH EP-ICIENCY CALCULATIONS 
NCOEF 4 NFTSt 10 t NSt 0 

c ,r: I Mrrn I aL DATA OF LOCAL CATCH EFFICIENCY I 


1 

-0.08 O.OA'5 


# *The names of the variables including those with $'s are printed 

to fill the appropriate spaces. Normally, these alphabetic symbols 
would not be input by the user but are given here to clarify the 
illustration. 



Figure 5.2 Particle trajectories released from x = -0.5, a = 0° 
(lower body) . 
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5.2 Axisymmetric Inlet at 30° Angle of Attack with Mach Number = 0.4 


The geometry for this case is the same as that shown in Figure 5.1. 
The input data files for generating the geometry with the program SCIRCL 
and axisymmetric potential flow program are the same as those given in 
Tables 5.1 and 5.2. The output data file Tape 07 of the potential flow 
program and Tape 05 (Table 5.5) created manually by the user are used 
as the input data files for the COMBYN computer program. The output 
data file Tape 17 (renamed as Tape 05) of the SCIRCL program, the output 
data file Tape 21 of the potential flow program, the output data file 
Tape 15 of the COMBYN program, and the data file Tape 02 created manually 
by the user (Section 3.2.3) are used as input data files for the trajec- 
tory program. As in Section 5.1, only four particle trajectories are 
computed. The trajectories are shown in Figure 5.3. 


TABLE 5.5 Input Data File Tape 05 for the Particle Trajectory Program 
(30° angle of attack). 


res'" CASE RELEASE 2.0 
2A2 90 242 90 1 0 i 0 

4-.4.00 14S.50 -60.000 ASS. ISO 0.00 0.000 

.000 .000 0.0000 .100 

2 81 90 -1 0 1 0 0 

.000 130.000 

0 . 9 P 3 2 0.2033 0 . 3 1 3 >’ 

.00 .00 0.50 S'J 


*The names of the variables including those with $'s are printed 
to fill the appropriate spaces. Normally, these alphabetic symbols 
would not be input by the user but are given here to clarify the 
illustration. 





6.0 DISCUSSION OF RESULTS 


This section discusses some of the limitations of the axisymmetric 
computer program in carrying out water droplet trajectory analyses. 

Section 6.1 describes modifications made to the original EOD program to 
accelerate the calculations of velocity at given points. Section 6.2 
describes a mesh generator program which was developed to investigate 
interpolation methods. Although the mesh generator when used in two- 
dimensional analytical solutions gave very good agreement, it was found 
to have no inherent advantages in terms of computational time for either 
the two-dimensional or axi symmetrical case. It is, therefore, not 
recommended and has not been incorporated into the trajectory program. 

This and other recommendations are described in Section 6.3. 

6.1 Modification to EOD Computer Program 

The original EOD program computed velocities at off-body points in 
parallel with on-body point calculations. All element contributions 
were calculated simultaneously for both on- and off-body points under 
control of subroutine MATRIX. The on-body point contributions were used 
to compute the system matrix which is then inverted while the off-body 
point calculations were carried forward for later computations. In 
subroutine PARTY each of the densities is read back from disk and all 
on- and off-body velocities are calculated in terms for the axisymmetric 
and crossflows. 

A straightforward application of EOD to calculate off-body velocities 
for the trajectories requires: 

1. Prespecifying all off-body points. 

2. Inversion of the system matrix. 

3. Multiple disk file reads to evaluate the off-body velocities. 

Repetition of Steps 1 to 3 is necessary to accommodate the trajectory 
solver which has no way to prespecify more than one point per trajectory 


before using that information to compute the next point, 

EOD was supplemented with a program which: 

1. Calculates the element contributions for one off-body 
point. 

2. Uses previous values obtained for densities thus 
eliminating matrix inversion. 

3. Retains all information required in memory during compu- 
tation thus eliminating a substantial portion of 
execution time needed for disk reads. 

EOD was modified to output the supplementary information needed by 
the trajectory solver. The supplement to the EOD code was then incorpo- 
rated in the trajectory solver. 

The calculation of velocity at a point is thus reduced to a straight 
forward integration. This integration is carried out by a Simpson rule 
algorithm which requires approximately one second* per point. Thus a 
typical particle trajectory calculation requiring approximately 500 
iterations to achieve the accuracy needed to compute local collection 
coefficients uses approximately 8 to 9 minutes of CPU time just to com- 
pute the velocities. Computing the particle trajectory between time 
steps requires at most 0.15 second (1.25 minutes per trajectory). Consid 
ering that roughly 10 to 15 trajectories are necessary to fully define 
the local collection efficiency, g, 2 to 3 hours of CPU time are used to 
determine one plot of g. Since the axisymmetric program consumes 85 
percent of the CPU time per calculation, a mesh generator with an inter- 
polation routine was developed and tested to determine if the economics 
of the computational time could be improved. 

6.2 Interpolation Method 

A sophisticated mesh generation program was developed in order to 
investigate the feasibility of using interpolation methods to find 
velocity components after computing their values at the nodal points. 

It was believed that greater efficiency could be achieved in the 


*A11 CPU times are based on a VAX 11/780 computer. 


calculation of trajectories if it were not necessary to calculate all 
velocities on a trajectory from the potential code (the direct method). 
The results, while promising, were mixed. 

For the two-dimensional flow case, the time required to generate 
the mesh was comparable to that required to run the specified test cases 
by the direct method. Once the mesh was generated, trajectory calcula- 
tions were approximately ten times faster than by the direct method. 
However, adding the time required to generate the mesh, the total compu- 
tation time of the local collection efficiency was essentially equal to 
and in some cases more than that of the direct method. 

Accuracy was thoroughly tested on a Joukowski airfoil using both 
the analytic solution and the potential flow solver (24Y). Relative and 
absolute accuracies in the velocity components of 0.1 percent were 
maintained throughout the field. 

For the axisymmetric case the mesh generation was considerably more 
time consuming. Even after eliminating matrix inversion, the potential 
code for this case (EOD) is a much slower procedure than 24Y consuming 
approximately one second per mesh point. Experience with the mesh 
generator indicated that approximately 17,000 points (17,000 seconds) 
would be required on the mesh to achieve 0.1 percent accuracy for the 
axisymmetric inlet case using the mesh generator. This gives more than 
five hours CPU time to generate the mesh. For the calculation of a 
local collection coefficient requiring a limited number of trajectories 
there is little question that the direct method is the faster procedure. 

6.3 Recommendations 

Based on the results of this study, the direct method of computing 
trajectories is recommended over the interpolation method. Both methods 
are extremely time consuming due to structure of the axisymmetric EOD 
computer program, despite the excessive modification to this program 
accelerates the velocity calculation. The EOD program is, therefore, 
not recommended as a viable computational tool for collection efficiency 
studies where numerous geometry configurations (whether due to design 
changes or due to ice build-up) are to be investigated. The physical 
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principles of the axisymmetric program are, of course, sound and fully 
tested. Therefore, if the program is to be used for icing droplet 
trajectory studies, it is recommended that the computational logic be 
completely restructured using modern developments in computer science. 
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APPENDIX 


PARTICLE TRAJECTORY COMPUTER PROGRAM LISTING 


DIMENSION Y< 6) t XD0T(6) t TITLE (20) »XLUOOO) » YL( 1000) t N( 10) » 
CFP( 10) »EP( 10) 


C-- PROGRAM FOR CALCULATING PARTICLE TRA JECTOR lES . 

rtKS UNITS ARE USED. 

C-- THE INITIAL INPUT OF PARTICLE DIA.-.ETER IS IN MIC?:ONS WHICH 

13 THEN CONVERTED TO «ET£R IN THE PROGRArt. 

C INPUT DATA FILES I 

INITIAL INPUT OF THIS FRAjcCTORY CODE 
AEROFOIL GEOMETRY. THE SANE FILE TAPE05 
USED IN THE NASA POTENTIAL FLOW CODE 

THE SAPlE AS THE OUTPUT F I LE ^ TAPE 2 1 OF NASA POTENTIAL 
FLGiJ CODE. USED FOP Air. vELuCITY CALCULATIONS. 

files : 

DATA STORED IN TPAEOl IS WPITTEN TO TArEOS LATTER. 
i'APE01 = TAPE03 FOR UNSYrt.'^ETRIC FLGU CASES 
DATA USED FOR THE CALCULATIG.S Or COLLECTION z.FFICI£nCY 
CREATED IF (LQPT.NE.O) . CONTAINS DETAILS OF THE 
TRAJECTORIES AS WELL AS INFORMATION OF WIND FIELD. 
CREATED IF(LQPT .NE. 0) 

DATA STORED IN TAPE08 IS FOR TRAJECTORY PLOTS. 

THE DATA STORED ARE XP.YP AND SW WITH THE 
F0RMAT(2X»3F10.S) . IF ( 3W . N£ . £8 . SBBS ) InPLIES 

THE END OF EACH TRAJECTORY. 

XP» YF' HAS BEEN NORMALIZED BY CHORD LEfJOTH. 

DATA FOR PLOTTING THE COLLECTIG.N EFFICIE.nCY 


COMMON/ COMB YN/VVREF 
COMMON/Q£iY/XA( 1000 ) » f A( 1000) 

COMMQN/XTN/XFRNT» YFRNT »XRERr YRERfNTCTALI 5) »NBDY . ID 
COMMON/XFR/XFCS) » Y F < 5 ) » X R ( 5 ) » YR ( 5 ) » Tnl K ( 5 ) . Cr MX < 5 ) f ( T MN •. 5 ) 
CQMMQN/SDS/S( 1 000 ) * SU ( 1 50 ) . ZO ( 1 50 ) » I E 

t.,DMM0N/IN2/XIN»YIN»£iFM»RL»VXR£r>VYRtF .PIT DOT tATK.r’IT .VREF » VXr'.v fF 

CQMMQN/F1/Q»CBAR» AMASS.Gf YYI . ALF A B » V I SCOS 

:aMMQN/MQD/LQPT»L£uM»XREAR»YLO»YUP 

COMMON/ I MM/ IP 1 » IP2 » NTQTM 

LCMMON/FRI/IFRNT (5) 

CCMM0N/DPD1/N3I » PL'JC (1 0 ) » I RS ( 1 0 ) » OPD < 1 0 ) .£PT 
CGMMQN/CDD/CF 

f(l)«X» Y(2)=Z» Y<3)*X1D0T. Y'4)sZlDC‘^» Y(5)*TH£T 

Y(6)=THETD 

REWIND 4 
REWIND 7 
R E W I N D S 

DATA PI/3 . 141592654/ 

R£AD( 2 » 252) N B D Y » NPL .• NSE AR . NEQ 
C“--“ CONTROL FLAGS 

R£AD(2»252) LEOM . LS YM . L R ANG 
CALL R E A D I N 

r;£AD( 2. 252 > L I M t LOPT . L PL C T . L TNL . LXOR 

C INITIAL CONDITIC.nS 

R£AD(2»254) O.VINF 
R £ A D < 2 r 2 5 4 ) D F' » P L * T I M S T ?' » Y L I M « C F 
READ ( 2 > 25 4 ) PRATK . PI T . PI TDCT » XORC f YORC 
READ<2»254) XREAR. YYLC. r r uP » YhA X . t M I r. 

R Eh £i (2.254) A[iEF3»ANLFL0.mNLCRD*EPSX 
READ(2.252) NSI 
IF(NSI .EQ. I )G0 TO 2 
DO 3 1*1. NSI 

F.£AD< 2. 221) PL WC ( I ) . DPD ( l ) ,CFP( I ) 

3 CONTINUE 
: bdyatk=anlflo 
ank=bdyatk+anlcrd 

=EI*PI/180. 

bOY ATK=BDYATK*PEI 

anlcrd=anlcrd»pei 

hNLFLO*ANLFLO<P£I 

ADE=-ADEFS 

,•1 D £ F' S = V I N F ♦ A Li E F' S 

vXREF^VINFtCGS ( ANLFLO > 

'J Y R E F * V I N F * S I N *. A N L r L 0 ) 

LIMA*LIM 

VREF = 5QRT (VXR£r*VXR£r+VYR£r*VYR'EF) 

;VREF*VREF 
VXP- VXREF 
VYP=VYREF 

.F(NBDY.£a.O)NBDY*l 

IF(ID.EQ.0)ID=1 

YnAX=YMAX«RL 

YMlN = YMlNi!RL 

ATK*PRhTN 

NSEA=NSEAR 

C SEARCH THE PROPER X-LOCATION TO RELEASE .-ARriCALS 

IF ( LXOR . EQ . 0 ) *30 TO 14 

NMAX*50 

XFRNT*XF( I ) 

XIN*XFRNT-0 . 2 
EPS=EPSX 

DY*< YYUP-YYL0)/10. 

NC=0 

134 CONTINUE 
NC*NCf 1 
DUX-O . 

/ IN*XIN-0 .2 

DO 131 1*1.11 

YP=YYL0 + FL0AT( I-l )tDY 


C 


TAPE02 

TAPE05 


TAPE21 


OUTPUT DATA 
TAPEOl : 


T A F' E 0 3 
TAPE04 


TAPEOS 


TAPE09 
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CALL VELCTY(XIN»YP»U1»U2) 
yT»30RT<Ul*Ul+y2»U2> 

D'J = ABS< I .0-VT> 

IF(DU.GT.DUX) PUX*D‘J 
131 CONTINUE 

IF(NC.GE.NHAX) GO TO 135 
IF(C:UX.GT.£FS) GO TO 13'i 
135 ’-ftIT£^6» 133) XIN»NC»NrtAX 
XIN=XIN«ftL 
GO TO 15 
•H XIN=XORC*Rl 
15 CONTINUE 

finish SEAFCHINo 

r-EI-FI/ieO. 
jF-rt = :.F*l .OE-OA 
AREA=0. 25*PI *DFn«DPM 
A H 0 F’ = I . C c 0 3 

ArtASS^PI *DPM«f3'«S:HOP/6 . 0 

VYI=0.4»AnAS3*<DPM/2.)**2 

CBAR=DPrt 

RH0=1 .23 

Q=0 .5*RH0»aREA 

^>'I3CQS-1 .5E~05 

yRITE<6*439)DPM»RL»TrMSTP»NPL»UXPEF.VYREF,AKK, PIT.PITD0T*G» 
lYYl f ADE 
NBY--NBDY 

IF'LTNL.Eu.O) GO TO 954 
f)BY=l 

ANLCRD=0.0 
954 CONTINUE 

IiO 94 IHK=1»NBY 
REWIND 3 
REWIND 9 
ID-IDK 

:f<ltnl .eg. 1) :d=i 

H=TIttSTP 
r PN T = XF t I D ) 

'FF:NTaYF< ID) 

.'iPERsXF; ( ID ; 

•F(XREAR .GT. XPEP) XREARsXPEP 
I P 1 = 1 

MTOTMsNTOTAL ' IDK) 

IFdOK .N£. 1) Ir l-HTCTAL ( IDK-1 > * 1 

:?2*IFRNT< IDK) 

rPEftaYRUD) 

’’HICKaTHIK' ID) 

YUP=YYrtX< ID ) +0 . 1 *THIK( ID) 

:L0*YYMN(ID)-0.i4THIK(ID) 

41j-A&3(ANLrL0) 
rPDUsABS < TP < : D ) -YF ( ! D > ) 

IF ;AQ.GT. 0.000001) THICKafrDWfiAR(ID)"XF<ID))*TAN<BDYATK) 
uRITE<6»27l ) 

.;RITE(4»272)NBDY»ID»XrRNT» YFRNT . <R£R » 'r PER , TH I CK » Y UP » YLO 
;F< ID.EQ . 1 ) 00 TO 9 7 
r.-!lNsZU«!RLt0.002 
YJIAXaYMINfO . 5tAB3 ( Yrt IN ) 

?7 CONTINUE 

DO 211 IJK=1»NSI 

LlM*LIrtA 

T3»0 

REWIND 1 

Ir^NSI.EQ.l) GO TO 4 
-RITE(6»437) UK 
DPh=DPD( IJK)*DP 
Cr=CFP( UK) 

DPM-DPM*! .OE-06 

AR_A=0.25»PI*DPH*DFrt 

AhASS = PI4iDPM*:«3*RH0P/6 . 0 

rYI=0.4*AflASS*4DPM/2. )**2 

CBAR=DPrt 

-.1 = 0 . S*RHO*AREh 

JRITE<6» 435) DPD( UK) f PLUC( UK ) »DPM»YYI 
4 CONTINUE 

IF<LEah ,£0.J)WRITE<6»136) 

WRITE(Af 433) 

COflPUTATION OF DROPLET TRAJECTORIES 

If LIM=1 SEARCH THE SURFACE RANGE WHERE THE PARTICLES IMF' ACT A 

IF<LRANG .NE . 0 ) CALL R ANGE ( YMAX r Yrt I N r T IMSTP * ADEPS r NEO ) 


YHIT=YMIN 
YMIS=YMAX 
IC = 0 

ITOl=100^ 

NSEAR=N3EA 
I rOT=NSEAR 
IPR0T=0 

IF(LI«. EO. I ) GO TO 41 

IT0T:=NPL*2 

IT01=NPL 

lUsNPL 

ZU=YYUP 

ZL= r YLO 

YUP=YYUP 

YLO=rYLU 

u iup=rup-i 


86 


a-a 


4 7 


as 


46 

48 



C 


12 


787 
798 
^ 75 


76 


79 


73 


C 

C 


214 


lUPsl SEARCHING UPPER LXHIT 
IUP»0 SEARCHING LOWER LIMIT 
CONTINUE 
TSP*T!.1STP 

Zr(LIrt .£0. 0) GO TO 47 
IF(SU( rU) .EQ.99.9999)IUaIU-l 
IF(ABSrt'.'IIS-YHIT) .LE. YLIM) GO TO 75 
YlN=<YHITfYMIS)/0. 

GO TO 48 

IF(NPL .NE. n GO TO 85 
YIN=rCRC*RL 
GO TO 46 
AC= I » 0 


DY=2U-2L 

Ir (LSYrt-.EQ. 1 ) DY = ZU 
DY-AC»DY /FLOAT (NPL-1 ) 
i'IN = AC«2U-Dr*<IU-iT0i) 

YIN=YIN*RL 

IF<IU .EQ. ITQT) GO TO 82 

IU-IU+1 

IPR0T=IPR0T+1 

IF( IPRQT.GT . 100) GO TO 743 
IFdU.EQ.l) YAD=YIN 
IF( lU.LE.NSEAR) GO TO 80 
URITE<6» 172) 

GO TO 236 

initially SU< I ) *88 .8888 

SU( I ) ^SURFACE DISTANCE IF PARTICLE HIT THE BODY 
SU(I)=99.9999 IF PARTICLE MOVES OUT OF RANGE 
SU(IU)*88.8888 

Y ( 1 ) = X I N 

Y < 2 ) - Y I N 

Y < 3 ) * V X P 

Y < 4 ) =t>YP 

r (5) =PIT»PEI 
Y(6)-PITD0T1!PEI 
:F(LEGM .n£ . I ) GO TO 92 

SET INITIAL CONDITIONS OF THE CASE FOR WHICH THE PARTICL 
15 EGUILI6ftl!JM IN THE AIR STREAM 
X*XIN/RL 
YpaY IN/RL 

CALL y£LCTY(X» YP»UX»UY) 

Y(5)*ATAN<UY/UX) 

Y(3)*UX*VR£F 

Y<4)=UYl:vREF 

Y<6)sO. 

CONTINUE 

ALFAB=ATK*PEI 

CALL ADAnS{Y»Xliur»TIhSTr'»ADEr‘S»IU»Yp*N£0) 

IF au .GT , ITOl .AND . SU< !U) .£0 .99.999?) GO TO 798 

iFCLIrt .£0. 0) gO TO 42 

SEARCHING SURFACE LIMITS 

IFrSWCU) .ES. J9.9999) GO TO 72 

YHITsY IN 

GO TO 42 

IFdU .EQ. 1) GO TO 73 

YHIS=Y In 

GO TO 42 

IW = 0 

IUP«2 

IF (YP.GT.YFRNT.OR.'rP.GT. fRER) GO TO 737 

rHIT=Y IN 

•30 TO 41 

CONTINUE 

YMIS*Y IN 

GO TO 41 


GO TO 236 

UPPER LIMIT GR LOWER LIMIT IS FOUND 
IC=ICf 1 

Ir<IC .EQ. 1) GO TO 76 
SL=SU( lU) 

ZL-=ZO( IW^ 

GO TO 78 
CONTINUE 
SU=SU( lU) 

ZU=ZO( lU) 

YMIS=tMIN 

YHIT*fAD 

IF(LSYrt.EQ.n GO TO 79 
GO TO 41 
WRITE- 6» 101 > 

SL-= -SU 
ZLf-ZU 

WRITE\6r321) ZUfSUfZLrSL 
IT01=IU 

E N El OF w I M I T '£ E A R C ri 

TOTAL COLLECTION EFFICIENCY 

IFCTHICK.EQ.O) GO TO 99 

£«<ZU-ZL)/THICR 

•FCNSI.GT. 1) EF(;JK)*e _ 

IF<NSI .GT. 1) URITE<6f 3r^)DPM»E 
IF ( NSI .GT . 1) GO TO 214 
WRITE(6.309) E 
CON r INUE 
WRIT£(6f 433) 


87 


32 

536 

575 

551 

552 


■vm 


234 


233 

211 


/ 2 1 3 
212 


1? 


C 

C 



IF(NPL.EQ.O) GO TO 236 
LIM-0 

itot=iu+npl 

nsear^itot 

2U*ZU-0. 00001 
2L*ZL+0. 00001 
GO TO 42 
CONTINUE 

IF<NPL.£Q.l) GO TO 19 

IA*0 

r3T = i 

ZMAX=Z0( 1ST) 

IflX^IST 
nsisT + 1 

liO 575 1 = 11 » ITOT 
ABZ = ABS<Z0( I ) ) 

IF (ABZ.LE. 0.0000001 ) Z0< I)=0.0 
IF<20( I ) .LE . ZMAX) GO TO S75 
IrtX^I 

ZMAX=Z0< I ) 

CONTINUE 
IF<IST .EQ 
IF(ZO( IMX) 

IF<SU( IMX) 

IF (LSYrt.NE , 1 .AND . Z0< I hX > 
IA=IA+1 

•JRITE< 1 ) 3U( IrtX ) f Z0< IttX) 
ZnUM=ZO( IhX) 

3DUM=SU( IMX ) 

Z0< IMX)=20< 1ST) 

3U(IMX)=SU(IST) 

Z0(IST)=2DUh 
3U< 1ST )=3DUn 
IST=rSTf 1 
IFdST .Lt. 

REWIND 1 
DO 232 I = lf lA 
READd ) SU( I) » Z0< I > 

CONTINUE 
.:s*iA 

Ir<LSYM.EQ.i) IS=IA» 

WRITE(3) IS 
DO 234 1=1 » lA 
URITE(3) SU<I)»20<I) 

CONTINUE 
•RS< lJK)»rS 

IF ( lA . £0 .IS) GO rO 211 


1) GO TO 551 

EQ . Z0< I3T-1 ) ) GO TO 552 
GE. 99.9999) GO TO 552 

.EQ.O.O) GO rO 552 


:T0T) GG to 586 


'1 


IA*IA-1 
DO 233 I=1»IA 
j = iA-n-i 
=UA=-SU( J) 

20A=-ZO(J) 

WRITE(3) SUAiZOA 
CONTINUE 
CONTINUE 
EPTsO.O 

IF<NSI .£0. i) GO TO 212 

DO 213 IK=lfNSI 

EPT = EPT + EP( iX ) *PLUC ( IK i 

CONTINUE 

CALCULATION OF LOCAL CATCH EFFICIENCY 
READ(2»120) ( TITLE( I ) . 1=1 »20) 

READ<2f252) NCOEF f NP TS f N S 

READ NASA/LEUIS DATA WHICH IS THEN STORED IN 

READ(2»120) ( TITLE< I) » 1=1 f20) 

r;EAIi(2»220) NASA 

■JRITE(9> NASA 

DO 17 I=lfNASA 

READ(2»22l)SF0IL»£FC 

WRITE<9) SFOILfEFC 

CONTINUE 

CALL £FFICY(NCOEF»NPTSf NS ) 

IF<LPLaT.EQ.O) GO TO 94 
IF(NPL .EO. 1) GO TO 18 
DATA PLOTING 

local CATCH EFFICIENCY • 

LDP = 0 
REWIND 9 
NA = 0 
NSET^^2 

DO 11 K=1»NSET 

READ(9) NX 

N(K)=NX 

IN=NAf 1 

NA=NA4NX 

DO 1 I=IN»NA 

READ (9) ;<L( I ) » YL ( I ) 

CONTINUE 

CALL APL0T(XL» YL »N»N3ET»LDP»FC»CZ) 

AIRFOIL geometry 

CONTINUE 

LDP = 0 


IFdDK .GT. 1) LDP = 1 
N( I )=NTOT 


CALL APLOT ( X A » Y A r N , 1 » L DP » FC » C2 ) 
TRAJECTORIES AND AIRFOIL GEONETRY 


TAPE09 


88 


CALL TRAJCT 
94 CONTINUE 
GO TO 236 
99 URITE(6»d6l) 

GO TO 236 
743 URITE(6»747) 

c formats 

101 F0RMAT(/»2X» 'SYMMETRIC FLQU CASE IS CCMFUTEH ' » / ) 

120 F0RMAT(20A4) 

123 F0RMATt/»2X» 'THE PARTICLES ARE RELEASEO FROM X/RL^ ' . 1F“E12 . 5 » / » 

I 2X»'UHICH 13 obtained AT THE'»I5»' LOOP OF'fiS.' LOOPS') 

172 FORMAT (/. 2X »' FORMAT 172:nO. OF PARTICLE TRAJECTORIES EXCEED 

i nsear» frogram stopped ') 

196 F0RMAT(/,2X» ' initial CONDITIONS FOR THE CASE L£0M=1'./) 

220 FORMATdS) 

221 FORMAT(3F10.5) 

252 FORMAT (S( 5X » I5»5X ). I 1 ) 

254 FORMAT(5(5X»F10.0) ) 

271 FORMAT (/ , 2X» ' N&DY' » ‘ ID ' . 12X» 'XFRNT ' » 10X» ' YFRNT ' » lOX. 'XREAR • . lOX » 

1 ' YREAR' » lOXf ' THICK' . lOX » 'YUPPP '• 10X» ' YLOOO' ) 

272 F0RMAT(2Xi2I4t6X»7(3X» 1FE12.5) 1/ ) 

309 F0RMAT(//f2X» 'TOTAL COLLECTION EFF IC lENCY ' » 3X » 1 PE 1 2 . 5 » / ) 

310 F0RMAT</»2X» 'TOTAL COLLECTION EFFICIENCY FOR THE DROPLET SIZE=' 

1 »E10.3» ' IS ' t£12.5»/) 

321 FQRMAT(/,6X» 'UPPER SURFACE L I M I T ' » 9X » ' LOWER SURFACE LlMIT'»/» 
l9XF'ZU'»12X»'SU'fl2X»'2L'*12X»'3L'»/»2X»4(E12.4»2X)»/»2X) 

4 33 FORMAT (SX »' XO ', ax » 'YO' f 23X »' XP ' »8X YP' »8X» 'SD' »3X »' DT ' ► 

1 10X» 'NSTP' ) 

427 FQRMAU /// 1 IX* ' C4«.F0LL0UING OUTPUT IS PQR DROPLET SHE 
1 DISTRIBUTION-- ',13*/) 

4 38 format </ » 2X»'DIAH£TER'»4X,' PERCENTAGE ' »4Xr 'PARTICLE' , 4X , ' INERTIA ' » 

1 /» 2X, 'DISTRIBUT. ' , 2X , ' OF W. CONT . ' * 3X * ' D I AME TER ' , 4 X , 

2 'MOMENT' ,//,2X,4( 1PE10.3,3X) »/) 

4 39 FORMAK / ,2X » 'PARTICLE' , 4X , 'CHORD' * 5X * ' T IME ' , 5X , ' MO . OF ' ,5X , 'UX' , 

1 3X» ■ OY 8X angle ' , 5X» ' PITCH' tSX, ' PIT ' ,7X, 'GRAVT '*5X, 

2 ' INERT. ',4X.'ADAMS',/*2X»' DIAMETER ' , 4X, 'LENGTH' * 4X, 

3 ' STEP' »5X, ' DROPLET ' , 4X, 'REF . ' ,6X» ' REF. ' ,5X» ' OF ATK . ' » 4X, 

4 ' ANGLE ', 5X, ' DOT 5X , 'CONSTANT ' » 4X* 'MOMENT' , 2X , ' CR I T ER I A ' , 

5 / / i3aPE? ,2 , IX > , 19 , IX ,3(IPE9 .2, IX) , / ) 

'47 F0RMAT(/, 'FORMAT 747: NO. OF TRAJECTORIES IS MORE THAN 

1 iOO, PROGRAM STOPPED ' > 

764 FQRMAT(/,2X» 'DATA STORED IN TAPE ' , I 3 » 2X , ' NO . OF DATA',I5>/ 

*6Xf 'S/RL' ,5X* 'Y07RL' ) 

793 F0RMAT(/,2X» 'Particle released between the region zl--zu 

: CROSSES UPPER OR LOWER IMPINGEMENT LIMIT TRAJECTORY, 

2 program STOPPED') 

361 F0RhAT(/,2X» ' THICK*0.0--PR0GRAM STOPED') 

:36 STOP 
END 

‘zUBRCUTINE ADAMS < ^L AST . YD 1 , TSP* EPS, I'J* YPL jNEQ) 


equartion solver 

THE AD AMS " nOUL TON r'RcD I C T OR -CORRECTOR METHOD IS USED 


DIMENSION YDl (6) , YLAST(6 ) , YFR£D(6 ) , TCCRR(6) fYP(4.6>»XPA(6),TPAd; 
DIMENSION A1 (6) ,A2(6) »A3(6) ,A4<6) 

COMMON/ XTN/XFRNT , yFRNT, XRER , Y RER • NTOTAL ( 5 ) »NPDY, ID 
C0MMQN/IN2/XIN,YINiDP,RL.VXR£F ,VYREF tPITDOT,ATKf PIT,UREF*UXP»VYP 
COMMON/ SDS/S( 1000 ) . 3U < 1 50 ) * ZC ( 1 50 ) * 1 3 
C OMMON/ MOD/ L OP T ,LE3M, XREAR ,YLO» YUP 

C IW,YPLr AND all THE COMMONS ARE USED EITHER 

C for program stop control OR FOR OPTIONAL PRINTS 

C LOOP form TIME3TEP TO START USING SIMPLE RUNGE-KUTTA METHOD 

H=rSP 
X=0 . 

IST0P*0 

NLOOP=lOO 

DX=TSP 

H024=DX/24.0 
DO 11 NT*1,4 
X 1 =x 

DO 1 J = 1»NECJ 
YPREDt J)=YLAST( J) 

1 YCORR( J)=YPRED(j) 

CALL F<Xl*YCaRR,YDl) 

DO 2 J-1, NEC 
Al<J)=H»YDl(J) 

2 YCDRR< J)*YPREDtJ)+AU J)/2. 

Xl=Xl+H/2. 

CALL F<X1,YC0RR,YD1) 

DO 3 J=L,NEQ 
A2< J)=H»YD1< J) 

3 YCORR< J)=YPRED( J)+A2< J)/2. 

CALL F(X1 , YCORR , YDl ) 

DO 4 J=1,NEG 
A3( J)=H»YD1 ( J) 

4 rCQRRC J) = tPR£D( J)+A3< J) 

Xl=XlfH/2. 

CALL r (XI fYCORRf YDi> 

DO 5 J=1»NEQ 
A4( J)=H4YDl ( J) 

5 YLAST ( J) =YLAST ( J ) + ( A I ( J J + 2 , 4 ( A2 1 J ) + A3< J ) 1 + A4 ( J i ) / 6 . 

A = X1 

DO 11 NVsWNEO 
11 YP(NT,NU)sfDl(NV) 

NSTPsO 
20 CONTINUE 


89 


31 

C 


42 


105 

?4 


jO 

100 


DO 50 IL00P*1 »NLOOP 

XsX+DX 

NSTPsHSTP+ I 

SIMPLE ADAMS PREDICTION 
DO 30 »NEQ 

YPRED<NV>=YLAST<NV)+H024*<55.0*YP(4tNV)-59.0<YP(3#Ny ) 
H-37.0*YP(2»N0)-9.0«YP(l»NO)) 

XPR=YPRED(1)/RL 

YPL=YPRED(2)/RL 

IrCXPR .LT. XFRNT) GO TO 111 

CALL MODE( XPRf YPL . ISTOP » lU. DX.NSTP) 

IFilSTOP .NE. 0) GO TO 52 
CALL F(Xf YPREDfYDl ) 

&ASFORD CORRECTION 

DLYMAX=0.0 

DO 31 NOsliNEQ 

YCORR(NO >sYLA3T(N0)+H024*(9.0*YD1(NV)+19.0*YP(4»NV) 
1-5.0*YP(3pNV)tYP<2»NO) ) 

IF(H0.NE.3.AND.N0»NE.4) GO TO 31 

CONTINUE 

IF<ABS<DLYMAX) .LE* EPS) GO TO 32 

HALVE DX 

DO 42 NV'IpNEQ 

DO 41 1=1 » 4 

XPA( I )=FLOAT( I ) 

YPA( I )=YP( I »NV ) 

YF<2»NV)=YP(3»NV) 

CALL TERP<2.5»YP1» YPPD»XPA» YPA»1»4»3»4) 

YP < 1 »NV ) =YP1 

CALL TERP<3.5»YP3»YPPD»XPA» YPA» 1 »4,3»4) 

YP(3iNy)=YP3 

CONTINUE 

NSTP=NSTP-1 

X=X-DX 

0X=O*50*DX 

H024=DX/24.0 

jO TO 50 

SHIFT VARAIAPLE FOR NEXT STEP 

CONTINUE 

DO 35 NV*1»NEQ 

r'P' 1 » NV) sYP(2 • NV) 

YP(2»NV)sYP(3»NV) 

YP(3fNV)sYP(4»NV) 

YP( 4pNV)*YD1 (NV) 

YLAST(NV)aYCORR(NV) 

PROGRAM STOP CONTROL AND OPTIONAL PRINT CONTROL 

XF *>LA3T< I )/RL 

YPL=YLAST(2)/RL 

CALL MODECXPtYPLr ISTOPt IU» OXfNSTP) 

:?•: ISTOP.NE.O) GO to 52 
IF(L0PT.£Q,0) go TO 94 
URITE(8»117) XP.YPL»SU(IU) 
r0RMAT(3F10.5) 

CALL VELCTY (XP.YPL»U1»U2) 

UX=-U1*VREF 

UZ=-U2*VREF 

<YLAST<3)+UX)«^2.0+< YLAST(4)+UZ)**2.0)*»0.5 
ALFA=ATAN<(-UZ-YLAST<4))/(UXfYLAST(3))) +YLAST<5) 

£1=YLAST(1)/RL 

£2=YLAST<2)/RL 

E2=TLAST(3)/yR£F 

E4=rLAST(4)/VREF 

£5=V/VKEF 

URIT£< 4 , 105) X»E1 »E2» YLAST( 5 ) » E 3 i E4 , YL AST ( 6 ) • E5 » U I r ‘J2 
F0RMAT(10E10.3) 

CONTINUE 

«»««««« END OF STOP control ********* 

GO TO 20 
CONTINUE 
JRITE<6f 100) 

STOP 

F0RMAT(6X» 'TIME-HALVING LOOPS EXCEED 100') 

RETURN 

END 

SUBROUTINE RANGE (YMX» Y MN * T I MSTP » ADEP5 f NEQ ) 

SEARCH THE UPSTREAM Y-RANGE SUCH THAT NONE OF THE PARTICLES 
OUTSIDE THE RANGE WILL HIT THE AIRFOIL. 

DIMENSION Y0(2) » YL (2) f Y< 6) fXDOT ( 6) 

C0rtnaN/IN2/XINfYIN»DF »RLfVXREFfVYR£FrPITOOT»ATKrPITfVREF»VXP» VYP 
COMMON /XTN/XFRNT » YFRNT » XRER » YRER » NTOT AL ( 5 ) tNBDY » ID 
COMMON/S DS/S ( 1000) »SU< 150) f Z0< 150) f IS 
C0MM0N/M0D/L0PT»L£C?M»XREAR f YLOrYUP 
PI=3. 141592654/130. 

YQ(1)=YMX 

Y0(2)=YMN 

IAC=0 

DY-l .0«( Y0< 1 )-Y0(2) ) 

IU = 0 
IPR0T=0 
I IST0P=0 
:w=rw+i 
IPR0T=IPR0T+1 
IF ( IPRQT.GT.30) GO TO 33 

su ( I w ) « 89 . sees 


L17 


90 


YIN=Y0( lU) 

r< 1 )sXIN 

Y(2)=YIN 

Y<3)«VXP 

Y< 4)=VYP 

Y(5)=PIT*PI 

Y(6)=PITD0T*PI 

IF(LEQrt.NE.l) GO TO 92 

C SET INITIAL CONDITIONS OF THE CASE FOR yHICH THE PARTICLE 

C IS EQUILIBRIUrt IN THE AIR STREAM 

XINN»XIN/RL 

YINN=YIN/RL 

CALL VELCTY(XINN»YINN»UXiUY) 

Y(5)=TAN(UY/UX) 

Y<3)=UX*VR£F 

Y(4)=UY»yREF 

Y(6)=0. 

92 CONTINUE 

YKL-Y0(2)/RL 
YKU=YO( 1 )/RL 
IF ( YKU.GT . YUP) YUP=YKU 
IF( YL0.6T.YKL) YLO=YKL 

CALL ADAMS (Y » XDQ T # T I hST P » ADEPS » lU » YP # NED ) 

YL< IU) = YP-YRER 

IF(SU< lU) .NE.99.9999) GO TO 10 
IF< lAC.EQ. 1 ) GO TO 29 
IF<IW.LT,2) GO TO 1 
29 YY=YL( l)*YL<2) 

IF< YY.LE.0.0) 60 TO 13 
IF<YL(1) .LT.0.0) GO TO 11 
GO TO 12 

10 IF(IU.EQ.2) GO TO 12 

11 Y0(1)=Y0(1)+DY 
IF(IU.EQ.2) IAC=l 
IU = 0 

GO TO 1 

:2 YO(2)=yO( 2)-DY 
IAC-0 
iUal 
GO TO 1 
13 YMXaYO(l) 

YMN»YQ(2) 
rO<l)*YO(l)/RL 
Y0<2)»Y0<2)/RL 
yRIT£(6»21) Y0(l).Y0(2) 

GO TO 34 

33 UR!T£<6»24) 

STOP 

. 21 FORMAT< /,2X» ' YMAXs ' . 1PE12 . 4* 3X» ' YMINs ' » IPE12. 4 ) 

24 FOfi;MAT(//»'IPROT .GE. 31 IN SUB^RANGE - PROGRAM STOPPED') 

34 RETURN 
END 

SUBROUTINE F(T»X.XDOT) 


C-- GOVERNING EQUATIONS OF THE PARTICLE MOTION 

^ — _ 

DIMENSION X(6) »XD0T<6) 

C0MM0N/IN2/XIN»YlN»DPfRLfVXREF»VYREF»PITDOTfATK»PIT»VREF»VXPrVYP 

C0MM0N/F1/Q»CBAR» AMASS»Gr YYI , ALFABfVISCOS 

COMMON/CDD/CF 

C XnOTd )=X1D0T» XD0T<2)-ZlD0Tf XDOT < 3 ) =X2D0T 

C XDOT(4)=Z2DOT» X DOT ( 5 > = T HETD f X DOT ( 6 ) = TH2D0T 

PI=3. 141592654 
XP=X< 1 )/RL 
i'P = X<2)/RL 

CALL VELCTY ( XP . YP » W 1 f U2 ) 

UX=-VREF*U1 

yZ«-VREF«U2 

VX=X(3)+UX 

VY-X(4)+UZ 

V=SQRT( VX«VX+VY»VY) 

RE = V1iDP/VISC0S 

CALL COEFF<RE»CDf CLfCM) 

CD=CD/CF 

IF(CL .£Q. 0.0 .AND. CM.EQ.0.0) GO TO 20 

IF(VX.EQ.O.) GO TO 17 

GAMMA-ATAN<VY/VX) 

IF<VX.LT.O.) GAMMA=GAMMA+PI 

IF( VX.GT.O.O.AND.vY.LT.0.0) GAMhA*GAMMA+2 . 4PI 
GO TO 19 
17 GAMMA=0.5*P1 

IF(UY.LT.O.) GAMMA=l.5*PI 
19 CONTINUE 

ALFA»X(5)-GAMMA 

ALFAD=ALFA-ALFAB 

alfab=alfa 

AL-Q»V4V«CL 
AD«Q*V*V»CD 
AM = Q*V*V4tCMtCBAR 
XDOT< 1 ) =X< 3) 

XD0T(2)=X(4) 

<D0T<5) =X(6) 

CG AMMA=CQ3 ( GAMMA ) 

SGAMMA*SIN( GAMMA ) 

XDOT { 3 ) =-( AL/AMASS ) * SG AMM A- \ AD/ AM ASS > <CGAMnA 
XDOT ( 4 ) s ( AL/AMASS ) tCGAMMA- ( AD/ AMASS ) *SGAMMA-G 
XD0T(6)=AM/YYI 
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GO TO 21 

20 CONTINUE 
XD0T(l)sX(3> 

XD0T(2)*X(4) 

XD0T<5>sX(6) 

XD0T{3)=-CD*C*U«VX/AMASS 

XD0T(4)*-CD*0*U«VY/AHASS 

XD0T(6)-0. 

21 RETURN 
END 

SUBROUTINE CUEF F < RE » CD » CL » CM ) 

C 

C-- COEFFICIENTS OF AERODYNAMIC FORCES 6ND OF MOMENTS 


C 

C 

c 


CL*0 . 

. 

IF(RE .EQ. 0.0) GO TO 22 
CD==24.0/RE 

IF(RE.LE.0.05) GO TO 14 

IF<RE.L£.3.0) GO TO 2 

IF(RE.6T.330) GO TO 5 

AQ=-28.339 

Al=38.969 

A2 = 0. 73204 

A3=-0. 00056084 

GO TO 10 

2 AQ = 0.0 
Ai=24 .167 
A2=3.2540 
A3 = -0 .23564 
GO TO 10 

3 AQ=0.0 
Al-93 . 462 
1^2*0,37576 
A3 = 0.0 

10 R2=RE*RE 

CD-<A0+A1«R£+A2«R2+A3*R2»RE)/R2 
JO TO 14 
22 CD5>0.0 
14 RETURN 
END 

SUBROUTINE rtOD£( XF »YP.ISTCP.IU»T»NSTr*) 


DECIDES UETHER THE PARTICLE HIT THE AIRFOIL OR NOT 


COMhON/SDS/3C 1 000 ) » SU < 1 SO ) » 20 ( 1 SO ) » 1 3 

CartMON/XFR/XX( S) »YY(S)»xB<5)tYB(S)fTHIK(5)fYYM:<(S>»YY?“H5) 
CGrtriON/IN2/XIN. rl N » DP f RL » VXREF » VYREF » r I TDOT f A TK t P I T . VREF f vX P , 
COMMON/OB Y/X ( 1 000 ) » Y ( 1 COO ) 

COrtMON/MQD/LOPT»LcQM»XReAR»YLOrYUP 

COMMON/ X TN/XFRNT »YFRNT » XRER , YRER » NTOT AL ( S ) »NBDY» ID 
COMMON/ IMM/ IP 1 » IP2»NT0TM 
URITE<8»«) » ' XP.YP.NSTP' »XP»YP.NSTP 
xFs(XP-XFRNT)*RL 

!F<A&S(XF) .L£. 0.00001) XPsXFRNT 

IC»0 

!FLAG=0 

::=o 

TF(XP.LE.XFRNT) GO TO 27 

•.-•:XP .GT. XREAR .OR. YP ,3T. YUF .OR. YP .LT. YLO) GO TO 32 
-- FIND X-POSITION OF PARTICLE RELATIVE TO X(I) AND X(IL) 
li IC-IC+1 

!r(IC .GT. 1) GO TO 7 
Ir(XOLD .GE. XFRNT) GO TO 5 
IFLAGsl 

YT-< YP-YPOLD)/ < XP-XOLD) * < XFRNT-XOLD > -f YPOLO 
IF<YT .GT. YFRNT) GO TO 3 
GO TO 6 

5 IF(YPOLD .GT. YFRNT) GO TO 3 
GO TO 6 

7 IF'YP .GT. YFRNT) GO TO 3 
A CONTINUE 

DO 2 IP*IP1»IP2 

I=*IP2-IP+IP1 

IL=I-1 

£F<IL .EG. 0) GO TO 27 
XR=< XP-X (!))*< XP-X( XL ) ) 

IF(XR .LE. 0. ) GO TO 29 

2 CONTINUE 
JO TO 27 

3 CONTINUE 

DO 4 IP*IP2»NT0TM 

IL=IP 

I-IP+1 

■:F(I .£Q. NTOTM+n GO TO 27 
XR=( XP-X( I ) ) *( XP-X( IL ) ) 

LF( :<R .LE. 0. ) GO TO 29 

4 CONTINUE 
27 YREF=YFRNT 

:F\YP .EQ. yfrnt.and.xp.eq.xfrnd go to 51 
rPT=(YP-YREF)/<rREF-YLO) 

GO TO 33 

29 YREF=<{Y(I)-Y(IL))/(X(I)-X(IL)))*<<P-<(IL))fY(IL) 

[F( IC .NE. 2) GO TO 31 
{I=ri+1 

{'EE=ABS(YR£F-YP) 

IF<YE£ .LE. 0.00001 ) GO TO. 51 


vYP 
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IF<11 .GT. 20) 60 TO 51 

XNEUsXP 

YREFN^YREF 

IFdFLAG .EQ» 1) GO TO 9 
GO TO 34 

- FIND THE PARTICULAR TRAJ. ACROSS THE BODY SURFACE OR NOT 
1 YPT*(YP-YR£F>/ (YREF-YLO) 

IFdFLAG .EQ. 0) GO TO 8 

IF(YT .GT, YFRNT .AND. YPT.GT. 0.0) GO TO 33 
IF(YT .LT. YFRNT .AND. YPT.LT. 0.0) GO TO 33 

3 XNEUsXP 
YNEUsYPT 
YREFN^YREF 
YRaYNEy*YOLD 

- FIND INTERSECTION POINT OF SPECIAL CASE 
IF( IFLAG .EQ. 0) GO TO 86 

9 A=<XNEU-XFRNT)*< YP-YPOLD)-<YREFN-YFRNT)*<XNEU-XOLrO 
B*YP-YPQLD 
C=YREFN-YFRNT 
D-XNEU-XFRNT 
E-XNEU-XOLD 

XPs(B»D»XOLD-C*E*XFRNT+< YFRNT-YP0LD)«ETD)/A 
YP=( D*B»YFRNT-£»C»YPOLDf <XOLD-XFRNT)*e»C) /A 
GO TO 1 

1 IF<YR .LE. 0.0) GO TO 34 
5 XQLD-XP 

fOLD^YPT 

rREFQsYREF 

ypOLDsYP 

ISTOP-0 

RETURN 

- FIND INTERSECTION POINT OF GENERAL CASE 

4 XP=<XNEU-XOLD) »( YPOLD-YREFQ ) /( YREFN- YREFO- YP + YPQLD > f XOLD 
YP=< YP-YPOLD) t < YPQLD-YREFO >/< YREFN-YREFO-YP+YPOLD) 4-YPOLD 

L CONTINUE 

IFdC.EQ.2) 1C=IC-1 
GO TO 13 
L YPaYREF 
GO TO 55 

- PRINTS DATA IF PARTICLE MOUES OUT OF RANGE 

2 XN=XIN/RL 
YN«YIN/RL 

yRITE<6»53) XN»YN#XP»YP>T»NSTP 

5 F0RMAT<2F10.5» ' OUT OF RANGE d 2F 1 0 . 5 » 1 OX . IPE 1 2 . 4 • I 1 0 > 

IST0P=1 

SUdU)s99.9999 
Z0( IU)aYN 
■30 TO 57 

- END OF DATA PRINT IF PARTICLE MOOES GUT OF RANGE 

- INTERPOLATION OF SURFACE DISTANCE 
i CONTINUE 

DO 24 IQaIPl»NTOTM 
1*10 
IM=I+1 

IF(YP .GT. YFRNT) CO TO 25 
GO TO 21 
5 I*NT0TM+IP1-IQ 
IM=I-1 

I XR1=(XP-X< I ) )«<XP-X< IM) ) 

IF<XRl.GT.O.) GO TO 24 

.-=(Sd)-3dM))/<X(I)-X(IM))«<XP-X<IM))+SdM) 

GO TO 59 
» CONTINUE 
? CONTINUE 
YN-YIN/RL 
XN=XIN/RL 
sudu)=p 

IS=IS+1 
Z0< IU)-YN 

URITE<6»56) XN»YN>XP»YPfPfT»NSTP 
S F0RMAT(2Fl0.5f ' HIT BODY AT ' r 3F 1 0 . 5 dPE 12 . 4 r 1 1 0 ) 

IST0P=1 
r CONTINUE 
RETURN 
END 

SUBROUTINE READIN 
DIMENSION A<500) fKT<2> 

COMMON/COMBYN/VOREF 

COMMON/COMBYNl/ JJr.NHUBPl f NHUBMX » AL IL » BL IL » RHOTOT r ATOT AL » 

1 X0N(500 ) f FACT0(2) f CCOSK 2) »CSINT(2) r 

2 Al(500»2) 

C0MM0N/XFR/XF(5) »YF(5)fXR(5)fYR<5)fTHIK(5) » YYMX( 5) f YYMN(5) 
C0MM0N/FRI/IFRNT<5) 

COMMON/OB Y/X< 1000) »Y< 1000) 

DIMENSION XX(SOO) »YY(500) 

COMMON/SDS/S< 1000) » SW < 1 50 ) » ZO (1 50 ) d S 
COMMON/XTN/XFRNT»rFRNT»XRER»YRER»NTOTAL(5)fNBDY,Z£iD 
COMMON /BLOCK!/ I 6E0MF ( 9 ) r I S I GF < 9 ) 

COMMON /BL0CK2/ HCURV ( SOO ) r HARC < 500 ) 

COMMON /BL0CK3/N0AXI » NOCROS t NOVORT » NOV I » N0V2 1 F IRSTE f L AST E » LSSE » SO * 
INSMALL 

COMMON /BL0CK4/A11 (500) » A12( 500) »A1 3(500) t A21 (500) » A22(500) » 

1A23( 500) 

C0MM0N/FL6/ HEDRdS) tCASE *NB »NNU 

lfFLG03 »FLG04 .FLG05 .FLG06 »FLG07 

l»FLG08 fFLGO? fFLGlO rFLGll »FLG12 

lfFLG13 #FL614 »FLG15 »FLG16 fFLG17 
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1»FLG18 fFLGl? #FLG20 tFLG21 *FLG22 

liFLG23 »FLG24 .FLG23 »FLG2A fFLG27 

COhnON/NNN/ NTi NO<U)» nN> NUNA(S)f 

INERlr NER2t NhA» NS!GA» NS1GC> 

1NUNC(5)» TYPEC(5)» NLF(H)» IEC» NSIGEC» 

1TYPEEC(5) »NUNEC<5) 

COMHON /KITZEL/ VX A < 5 ) * UY A < 3 ) * SIGA(500t5)» 

1 yXC<5>fyYC<5) .UZC(5) »3IGC<500»5> » 

2 yN(5).VT(5) 

DOUBLE PRECISION HEDRt CASE 


TYPEA{ 5) » 


FLG03 

»FLG09 

»FL6i4 

»FLG19 

,PLG24 


»FLG04 

tFLGlO 

.FLG15 

»FLG20 

iFLG25 


»rLG05 
•FLGl I 
»FLG1S 
.FLG21 
»FLG26 


.FLG06 
»FLG12 
fFLGl? 
f FLG22 
•FLG27 


f FLG07 


X2(500)f Y2(500)» DELS(500). 


INTEGER 
IfFLGOa 
1 f FLG13 
1 f FLG18 
l*FLG23 

REAL MNiLSSEfMX.MY 

INTEGER BDNiSUPKS 

COhMON /CL/ Xl(500)f Yl<500)» 

1SINA<500) fCQSA<500) 

LOGICAL NOAXI f NOCROSf NQUQRT f NOVI fNQV2»FIR6TE»LASTE 

C READ INFQRHATION FOR CALCULATING WIND VELOCITIES 

EPS=0. 0000001 
DR=0. 0174532925 
REUIND 15 
READ< 15)NTHETA 

READ(15) JJfNHUBPlfNHUBttXf ALlLfBULfRHOTOTf ATOTALfTHEAfXON 

DO 12 I=1»NTH£TA 

R£AD( 15)FACT0RfCCQSTHf CSINTHf A 

FACTO< I >sFACTOR 

CCQSTCD-CCOSTH 

CSINT(I)=CSINTH 

DO 11 J=l»500 

11 Al( J,n=A< J) 

12 CONTINUE 
REUIND 21 

R£HD<2UIG£QrtFf ISIGF 

READ<21 )X1 > Y1 f X2f Y2f DELSf SINAfCOSA 

READ ( 21) HEDRf CASE fNBfNNUf 

1 FLG03f FLG04iFLG05tFLG06f FLG07»FL608f FLG09f FLGlOf FLGU r 

2 FLG12»FLGl3f FLQ14»FLGl5»FLGl6»FLG17»FLGl8fFLGl9fFLG20f 

3 FLG21 »FL622f FLG23f FLG24.FL625fFL625fFL627 
READ(2l)NT»ND»NNfNUNAf TYPEAfNERl » NER2 f NrtA f NS I GA » NS I GC t NUNC f 

1 TYPEC^NLFf lECfNSIGECfTYPEECfNUNEC 

READ(2UHCURVfHARC 

R£AD<21 )NOAX I .NOCROSfNOVORTf NOVI fN0V2. FIRSTS f LASTS f 
i LSSEfSOfNSMALL 

REAIK 21) A1 1 f Al2f A13 f A21 f A22» A23 
READ<2i)SIGAfSIGC 
TYPEFf 'FINISH READING ' 

REUIND S 

C READ X-f Y-i COORDINATES 

C READ TITLE 

NTOTsO 

READ<S»SO)HEIiRf CASEf PSFf NB 
50 F0RMAT(lSA4f5X.A4f7XfA4f/f28Il) 

READ<5»60)CH0RD»HNf TCNSTfEPSLON 
60 F0RMAT(4F10.0) 

Nll=l 

N12*0 

DO 30 K-lfNB 

SEAD(5f 10) IGEOMF<K ) f ISIGF<K ) f I CURVN » NONEUF f IFORHT fNN.HX . 

I MYf THETAf ADDXf ADDY 

iO F0RMAT(5I1 » ISf 5F10.0) 

READ(5f 20) BDNf 3UBKS »NLF(K ) f XE» YE 
20 F0RMAT<3(5Xf I5)f2F10.0) 

N22=N12+NM 

IF( IF0RMT.EQ.2) GO TO 120 
:F( IFORMT»EQ. I ) GO TO liO 
C IF0RMT = 0 

READ<5f710)(XX(I)fI-Nll»N22) 

READ<S»710) < YY ( I ) r I=Nll f N22> 

GO TO 130 

110 READ(5»790) (XX<I)fYY(I)fI=Nll f N22) 

GO TO 130 

120 READ(5f800) <XX<I)fYY(I)fIsNl 1 f N22) 

710 F0RMAT(6£13.4) 

790 FORHAT<2F10 .5) 

300 F0RHAT<F10.5f lOXfFlO.5) 

130 CONTINUE 
Nll=Nll+NN 
N12=N12+NN 
KT<K)=NN 
30 CONTINUE 

NTOTAL( I )=KT( 2) 

NT0TAL(2)=KT(2)+KT(1)*2-1 
NT0TAL<3) =KTt 2) +NT0TAL ( 2 ) 

DO 35 KalfN22 
X(K)=XX(N22-K+l) 

Y(K)s-YY(N22-K+l ) 

IK=N22fK 

IF(K.EQ.N22) GO TO 35 
X< IK)=XX<K + t ) 

Y( IK >«YY (K-H ) 

35 CONTINUE 

DO 41 K=lfNBDY 
IIB=NT0TVI 
I IA*NT0T+2 
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NT0T=NT0TAL(K) 

3(IIB)*0.0 

<XO-X(II&) 

NN0*1 

DO 2 I=IIA»NTOT 
DDX*X< I)-X ( I-l > 

DDY=Y(I)-f(I-l) 

S ( I ) sS ( I- 1 ) +SORT ( DDXf PDX+DDY*DDY ) 

IP( X( I ) .GT. XXQ) GO TO 2 
XXO=X< I ) 

NNO=I 

2 continue 

SC*S(NNO) 

DO 4 I=IIB»NT0T 
3(I)=30-S(I) 

4 CONTINUE 

IF( A&S<rtX ) .LE . EPS) GO TO 200 
DO 210 I=IIB»NTQT 
210 X<I)=X(I)*rtX 

200 IF(ABS<MY> ,LE. EPS) GO TO 220 
DO 230 I=IIB»NTQT 
230 Y(I)=Y(I)»MY 

220 IF< ABS<ADDX) .LE.EPS) GO TO 41 
DO 270 I=IIB.NTOT 
270 X(I)=X(I)+ADDX 
41 CONTINUE 

CALCULATION OF SURFACE DISTANCE 

lAsl 

NTOT-NTQTAH 1 ) 

XMAX-X( 1) 

XOO=X(1) 

YMIN=Y( 1 ) 

YhAX=Y(l) 

DO 91 K=*1»NBDY 
IF( K .EQ . 1) GO TO 51 
IAsNT0TAL<K-1)’K 
NTOT*NTOTAL(.\> 

<flAX = X< lA) 

;<00*X( lA) 

YNAX=Y< lA) 
rMIN=Y<IA) 

31 IBsIA+l 
N0»1 

DO 21 lalBfNTCT 

:F(X( I) .GE.XOO) GO ro 3l 

:<oo=x < I ) 

No*r 

GO TO 21 

SI IF(X< I) .LT.X.1AX) GO TO 2l 
XNAX*X< I ) 

NRsI 

21 CONTINUE 
IFRNT<K)aNO 
XF(K)*X(NO) 

YF(K)sT (NO) 

XR<K)*X(NR) 

rR(K)=Y(Nft) 

DO 61 I-IBfNTOT 

IF ( Y C I ) .GT . YrtIN) GO TO 62 

YMIN=Y< I ) 

GO TO 61 

62 IF(Y( I) .LT. YNAX) GO TO 61 
YMAX = Y( I) 

61 CONTINUE 

THIN(K)=YMAX-Y«IN 
YYMX(K)=YMAX 
rYMN(K)=Y«IN 
51 CONTINUE 
RETURN 
END 

SUBROUTINE UELCTY<XP, YP f UXCC » U YCC > 

COHMON/COMBYN/UUREF 

COnrtON/COrtSYNl/ JJ»NH'JBP1 f NHUBMX» ALIL ♦ BL IL » RHOTOT » A TO T AL » 

1 X0N(500) f FACTORS 2) •CCDSTH( 2) *CSINTH( 2) • 

2 A(500f2) 

COMMON /KIT2EL/ VXA < 5 ) f U Y A ( 5 ) f SIGA(500.5)» 

1 UXC(S)fUYC(5)»VZC(5) fSIGCCSOOfS) f 

2 VN(5)fUT(5> 

UREF=VVREF*100. /30. 48 
XPKT-XP 

Y PM=ABS(YP) 

CALCULATE THE INDIVIDUAL SOLUTIONS 

CALL MATRIX(XPKTf YPKT) 

CALL AXIS 
CALL CROSS 

TO CALCULATE THE COMBINATION SOLUTION 

IT = 1 

IF(YP.LT.O.O) IT=2 

EXCON=l .0954«ATOTAL»125./216. 

V30NCCsAT0TAL/3QRT( 1.2) 

IF(XPKT.LE.X0N( JJ) )G0 TO 115 

CALL S INTP( XON ( rmUBPl ) f A< NHUBP 1 * I T ) f J J-NHUBMX f XPKT f AR > 

(30 TO 120 
115 ARsA ( JJ f-I-T ) 

120 OB^FACTORC IT)/AR 

VXCC=ALIL«VXA( I )+PLIL*VXA( 3) +CC0STH( IT) «UXC( 1 ) 

V yCC = AL IL«Ur A ( 1 ) f SL IL*VY A( 3 > +CC OS TH( I T) *VYC C I ) 
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UZCC»CSINTH< IT) tVZC< 1 ) 

CALL VBARIT( AFOTALf RHOTOr iRHB) 
yRES0F = 3£)RT ( yXCC»< 2 + V YCC*« 2 + yZCC » #2) 
R£SVV=VR£SOF 
RBFORTsRHB/RHOTOT 

*JR£SOF3yRESOF<( 1 . /RBFQR T ) ** < VRE SO? /‘J& ) 

yXCC*yRESOF»VXCC/RESW 

>;YCC=VRESOF*'JYCC/RESyV 

';ZCC = yRESOF*y2CC/RESVV 

IF(VRESOF.LT.VSONCC)GO TO 200 

VSAUE=ORESQF 

■;CONA= .2»< VSAOE/ATOTAL ) t»2 
PSOPTC=( 1 . -VCONA) *»3 . 5 
RHORTCsVSAVEcFSOPTC** . 71S/EXCQN 
rF(RHORTC.EQ.O. ) RHQRTC = 1 . 

VRESOF =OS0NCC* <l. + (ySAOE/’-'SGNCC-l.>*»<l 
CHANGE=URES0F/USAOE 
VXCC-OXCC^CHANGE 
OYCC=yyCC*CHANGE 
VZCCs02CC*CHANGE 
200 CONTINUE 

UXCC=yxCC/OREF 

uycc^oycc/uref 

IrvIT.EQ.2)VYCC=-VYCC 

RETURN 

END 


SUBROUTINE S I N TP ( Z » W » N » X 1 » Y 1 ) 
DIMENSION X(400)»Y(400)»Z(l>»U(l) 
DO 10 1 = 1, N 
X(I)=Z(I) 

10 Y(I)*W<I) 

CALL SORTXY (X,Y,N) 

DO 15 I= 1 ,N 


/RHORTC) ) 


K = I 

IF (Xl.GT.X(D) GO TO 15 
IF (Xl.EQ.X(I) ) GO TO 20 
IF <X1 .LT.X< I) ) GO TO 25 
15 CONTINUE 
20 Y1=Y(K) 

■30 TO 30 

25 IF <K*£Q.l) GO TO 35 
IF (K.EQ.N) K=N-1 
IF <X(K) ,£Q.X<K+1 ) ) K«K-1 

ai=<Xl-X(K) ) *(X1-X<K+1 ) )/(X(K-l)-X(K) )/(X(K-l)-X(K+l ) ) 
li2s<Xl-X<K-l))»(Xl-X<K + l))/<X(K)-X(K-l))/(X(K)-X(K + l)) 
^3*(X1-X(K-1))*<X1-X<K))/<X<K+1)-X(K-1))/<X(K+1)-X(K)> 
Yl=Y(K-l)»UlfY<K)»U2fY(K+l)*U3 
JO RETURN 
35 Yl»0.0 
aETURN 
END 


SUBROUTINE SORTXY ( X , Y , NPTS > 

DIMENSION X< 100) » Y( 100) 

10 N*NPTS 
15 NN*N-1 
20 DO 55 KT=1»NN 
XMIN=X (KT ) 
jaD=KT 
JKL=KT+1 

25 DO 45 JK=JKL»N 

30 IF (XMIN-X(JK)) 45,45,35 

35 XMIN=X^JK) 

40 JAD=JK 
45 CONTINUE 
50 YMIN=Y(JAD> 

A< JAD)=X(KT) 

Y( JAD)=Y<KT) 

X(KT)aXMIN 
Y<KT)=YMIN 
55 CONTINUE 
RETURN 
END 

subroutine UBARIT< UBAR, ATOTAL ,RH0T0T *RH0BAR) 


C APPROACH 5 

C TO SOLVE VBAR COMP ITERATIVELY 

C 

VCRIT=AT0TAL/SQRT( 1 .2) 

IrrO 

VGUES-sVBAR 

10 VGUESA=<VGUES/AT0TAD**2 
A=1 .0' .24VGUESA 
B=A-UGUESA 

OCOhP*^ VBAR-A412 .SiVGUES)/ ( A»«l . 5*B)+VGUES 
IF (ABS( <VCOMP-VGUES)/VCOrtP> ,LT. .0001 ) GO TO 15 
I = I-fl 

IF (VCOMP.GE.VCRIT) VCCMP= . 5< < VGUES-fVCRI T ) 

VGUES^UCOMP 

IF ( I .GT . 20) GO TO 15 

GO TO 10 

15 RHOBAR=( 1 .0- . 2* < VCOMP/ ATOT AL ) « «2 ) ▼»2 . 5 *RHO TO T 
IF (I.GT.20) URITE < A , 20 ) VB AR , VCOMP • PHuBAR 
IF(I .GT .20) V&AR=VCOMP*RHO&AR/RHOTOT 
20 FORMAT (1H0,34HI EXCEEDS 20 ITERATIONS FOR RHOB AR , 5X , 7H VB AR » »IPE 
il0.3,2X,8HVC0MP = » 1 PE 1 0 . 3 , 2X » 9HRH0BAR = -IPEIO.3/ 

270H VBaR has been REDUCED TO VCOMP*RHOBAR/RHOTOT , UHERE VCO«P=VCRIT 
ilCAL ) 
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RETURN 

END 

C*»»$#<*»«»a£GINNING Or ELEMENT ASINE/ 
function ARSINIX) 

ARSIN«AS1N(X) 

RETURN 

END 

C*»»*»»*«»»BEGINNING OF ELEMENT AXIS/ 
SUBROUTINE AXIS 

C t COMPUTE AXISYMMETRIC 

COMMON/NNN/ NTi ND(ll). 

iNEKl» NER2» NMA» 

iNUNC(5>. TYPEC<5)i NLFdl). 
1TYPEEC(5) »NUNEC<5) 

REAL MN 

COMMON /KIT2EL/ UX A ( 5 ) i U Y A < 5 ) t 






VELOCITY COMPONENTS AffD PRINT 
rtNf NUNAiS)* TYP£A(5>* 

NSI6A# NSIGC» 

I EC f NS I GEC t 


£IGA(500*5) • 


I 


:OMMON /TL/ 


vxC(5>»VYC(5)fVZC(5) »SIGC(500»5) » 
VN(5) »VT<5) 


A(500) I 
CX(SOO) f 


DX» 
XK » 


B( 500) * 
CY<500) • 
J » 

DY» 
cEK f 


L = 0 

DO 050 N>1»NSIGA 

SX»0 . 0 

SY= 0*0 

* NO, OF ELEMENTS LOOP 
DO 510 J=liNT 
SXsSX+AX( J)#SI6A( J>N) 

510 SY=SYVaY( J)«SIGA< J»N) 

IF ( N ,NE . 1 ) GO TO 520 
VXA(N) = 3XF1, 

GO TO 540 

320 IF (NUNA(N-1 ) ,NE. 

L*L+1 

vXA(Ni*SX+VN(L) 

VYA(N)=SY+VT<L) 

SO TO 550 
530 VXA^N) s SX 


AX(500) » 
C2(500) » 
J1 r 
NI f 
EKK » 


AY < 500 ) » AZ ( 500 ) » 
AXV(500) »AYV< 500) i 


S J > 
XJf 
K> 


DS » 
YJ» 
r'F 


.123456) 60 TO 530 


540 VYA(N) = 3Y^ 
CONTINUE • P#/ 

— n A 'A 


550 


GO TO 580 

■ ADJUSTMENT 


530 

5S0 


(MN.EQ.070) 

* nACH NO 
DO 560 N«1»nS1CA 
VYA(N)sVYA<N)/SGRT( 1 ,0-MNtfifO 
VXA( N ) s( VXA(N) - i . ) / < 1 .0-rN«r.N 
CONTINUE 
‘sETURN 
£ND 

l^*$»4**»«*|(£3inning of element CRQS/ 

SUBROUTINE CROSS 

C * COMPUTE CROSS FLOW 

■:ommon/nnn/ nt» ND<ii)f 

INERlf N£R2» NMAp 

lNUNC(5)f TYPEC(S)» NLF<U)» 
iTYP£EC(5) »NUNEC(5 j 


FI . 


VELOCITY 

NSIGAf 

I£C» 


components and print 

NUNAiS) • TYPEA<5) . 
NSIGCp 
NSIGEC/ 


REAL 
COMMON 


MN 

/MT2EL/ 


COMMON /TL/ 


VXA\S)»VYm(S). 
UXC(5) »VYC(3) pVZCd 
N<5) rVT<5/ 


SIuA(500»5) • 

) f 5IGC(500»5) . 


‘^OGICAL FF 
DO 370 N=1fNSIGC 
SX»0 . 0 
SY»0.0 
ipsO.O 


A(500) p 
CX(500) I 


DX P 
XK p 


B(500) p 
CY (500) • 
Jp 
DY p 
EEKp 


AX ( 500 ) p 
CZ(50O) 1 
J1 p 
NI p 

C.NK p 


AY(500) p AZ'SOO) p 
AXV(500) • AYV( 500) p 


Sjp 

xjp 

K, 


LSS P 

Y Jp 

DC- 


OF ELEMENTS LOOP 


« NO . 

DO 350 J=lfNr 

SX^SX+CX<J)»SIGC(JpN) 

5 Y 3 SY+CY( J)»SIGC( JpN) 

350 SP*SP+CZ( J)*SIGC( JpN) 

VXC(N)=SX 
VYC(N)=SY+1 . 
vZC(N)=SP+l . 

370 CONTINUE 
RETURN 
£ND 

C»:«»*«*»t4!»BE0INNING OF ELEMENT ELINT/ 

SUBROUTINE EL I NT 3 < XKSQ p XN p PH I p P I E > 

C THIS SUBROUTINE CALCULAFES THE INCO.MFLETE ELLIPTIC 
C THIRD KIND. THE ARGUMENTS ARE' 

C XKSQ value OF N SQUARED 

C XN VALUE OF MINUS ALPHA SQUARED 

C PHI VALUE OF PUl 

C PIE VALUE OF INCOMPLETE ELLIPTIC INTEGRAL OF THIRD KIND 

DATA HP /I .570796/ 

DATA ROUND /. 0000050/ 

3K=XKSQ 

FNsXN 

r=PHI 

•F (FN .EQ. -I . 0. ANP.sk .£Q. 1 .0) GO TO 480 
iF(SK*GT.l.) GO Tu 470 


iNrEGRAL UF THE 


ORIGINAL PAGE IS 
OF POOR QUALITY 
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c 

10 . 


20 

20 


40 


50 

60 

70 

ao 

?o 

100 


no 


120 


130 


c 


140 

130 

160 

170 

130 

190 


200 


IF(FN.LT. ( -1 . ) ) GO TO 470 
IF(P) 10 » 470 » 20 
NORMALIZE PHI 
A = -l , 
p = -p 

GOTO 30 
A*l . 

B-l . 

SB-l . 

IF (ABS<P-1 .570796 ) . LE . 1 0 . 0 * c ( - 7 > ) GO 

IF(P-HP) 110 p 100 » 40 

J-P/(2.*HP) 

,rx = 2* J 
P1-P-XX*HP 
P = HP 
B = ” 1 » 

GOTO 100 
0 = 3Urt 
8 = 0 . 

IF<P1-HP) 60 I 70 » 80 
P = Pl 
XXX*1 . 

GOTO no 
PIE=(XX+1. )*A»D 
GOTO 460 
XXX=-1 . 

XX=XX+2 . 

P=2.4HP-P1 
GOTO 110 

PIE=A»(XX*D+XXX*SUrt) 

GOTO 460 

IF(SK.£Q.l.) GOTO 470 

IF( FN . EQ. ( -1 . ) ) GO TO 470 

IF(P.GT. 10.E-4)G0T0 130 

IF<FN.GT.O. )GQTQ 120 

5Urt = P 

GOTO 440 

RRT=SQRT(FN) 

SUf1=ATAN<P*RRT) /RRT 

GOTO 440 

S=SIN<P> 

32=S*«2 

C=COS(P) 

:F<3K.GT.0.64)G0T0 210 
IF<ABS<FN> .GE.0.6)G0TC 160 
POUER SERIES IN N AND K SQUARED 
SA«1 . 
i&s3K/2 . 

CB»S*C 
C A * P' 

FMsO . 

£ u rt = ?■ 

.X»SUfi*l .E-8 
SA-SB-SA4FN 

C As(-CB/(2.*(FM+1. )))+(!. -.5/(FH+l.»)*C A 

t =3A4CA 

3UM=SUn+Y 

IF( CS&-4CA) .GT.X) go to ISO 
rF<ABS<Y) .LT.X) go TO 440 
Fh*FM+l . 

GB=CB*S2 

Ib-=(1.“T5/(F« + 1.))*SR«SB 

GOTO 140 

r'OUER SERIES IN K SQUARED 


ST-SGRT< l.+FN) 

IF<RT.NE.O.) GO TO 170 

G-S/C 

GOTO 190 

IF(C.GT.4,£-3)G0T0 180 
G=<HP-(C/(RT*S) ) )/RT 
GOTO 190 

G=ATAN<RT»S/C)/RT 
C1=G»1 .E-3 
■£ = P 
F = 3«C 
H=l . 

£Uh = 6 
FM = 0. 

G-<EtG)/FN 

H = H1f< 1 . -0.5/<FMf I . ) ) 

G2 = H4tG*Pf'. 

SUf1 = SUM + G2 

IF(G2.L£.G1)G0T0 440 
FH=Frtf 1 . 

£--F/(2.4FM) + U.-0.5/Frt;4E 
F=F*32 


TO 100 


c 


GOTO 200 
JIO SKP«1.-SK 

IF(S.LT.C) GO TO 320 
ADDITION FORMULA 
ZP=SORT' 1 .-SK*S2) 

RTl=3QRT(AbS(FN«<FN+l.)*<FN + Sh. )>) 

£3Ta(l.-ZP)/(SK<(C+l.)) 

tps ( SSTtS»RT I ) / < I . +FNtS2-FN «S3 r*C*ZP > 

iF(FN) 220 » 290 . 230 
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220 IF<RTl.NE.O.) GO TO 240 
RaS/(C+l.) 

IF(FN,NE. < -1 . ) ) GO TO 230 
CFa(2.«R-SK*SST*S-<S/C)*ZP)/SKP 
GOTO 300 

230 CF=c SST*(S“(2./R) )+S*C/ZP)*(SK/SKF) 

GOTO 300 

240 IF(FN*<FN+SK) .LT.O.) GO TO 260 
250 CF=(FN/RT1)»ATAN(XP) 

GOTO 300 

260 IF( ABS(XP) .GE.O . 1 )GOTO 270 
YX=XP#«2 

YX=2.*XF»<l.+YX»<l./3.+YX*(.2+YX/7.))) 

GOTO 280 

270 YX*ALOG< ( 1 .+XF)/( 1 .-XP ) ) 

290 CF=(FN*YX)/<2.*RT1) 

GOTO 300 
290 CF=0* 

300 

S=SORT<SST) 

C=SQRT( 1 .-SST) 

GOTO 320 

310 SUM-2. 4SUrt“CF 
GOTO 440 
320 U=3/C 
y=i ./C 
T=U»U 
U»U««2 

IF<S .GT .0. 1 )GQTQ 330 
R=S*(l.fS2*<l./3.tS2*(.2+S2/7.))) 

GOTO 340 
330 R = ALOG(U + »J) 

340 D=1.+FN 

IF(Ii.GT.SKP) GOTO 360 

C POUER SERIES IN IfN AND 1 - <K SQUARED) 

CAal . 

:Ba-O.S*SKP 
hL=(T+R)/2.0 
BE=U*U**3 
FM*0 . 

5UM=AL 

Tl*SUM*l .E-8 
350 CAa-D*CA+C& 

AL=(B£-< 2 . »FM+1 . )»AL)/ ( 2.<<FM + 2. ) ) 

X=CA*AL 

SUMaSUH+X 

IF ( ABS<X) .LT.TDGOTO 430 

FMapM+l , 0 

CBa-U2.*FM+l.)/(2.4<FM+l.)) )*CB»SKP 
IF <A8S< BE) .LT . 10. E- 30) B£aO . 

BE-&E«U 
GOTO 350 

C F'OUER SERIES IN 1 - (K SQUARED) 

360 RTaSQRT< ABS<FN) ) 

:F(FN) 370 * 380 » 390 
370 QaALOG< < 1.+RT«S)/<1 .-RT*S) )/(2.*RT) 

GOTO 400 
330 Q=S 

GOTO 400 

390 QaArAN(RT*S)/RT 
400 SUMaFN*Q+R 
PKpaSKP 
Hpa-0.5 
FMal . 

T1=SUM»1 .E-8 
410 a=(ft-Q)/D 

RaT/< 2. *FM)-( 1 .- .5/FM) *R 
XsAP»(FN»Q+R)*PKP 
SUMaSUMf X 

IF(ABS(X) .LT.TDGOTO 420 
T = T*U 

PKPaPKPf SKP 
FM = FM + 1 . 

AP=-AP*(1 .-.5/FM) 

GOTO 410 
‘120 SUMaSUM/D 
430 IF<BB.LT.O. )GOTO 310 
440 IF(B) 50 » 90 » 450 
450 PIE=A»SUM 
460 PIE=PIE+PIE*ROUND 
RETURN 

C ERROR RETURN 

470 PIEaO. 

GOTO 460 

2 CASE OF PI< D DPMI ) 

400 PIE a 0.5*<TAN<P )/ COS<P ) + AL OG < T AN ( ( HP+P 
GO TO 460 
END 

C*4!»*«j[$$«4t&EniNNING OF ELEMENT ELIFV 
SUBROUTINE EL IP ( XK » EEK » EKK ) 

C < START 

ETA s 1. - XK 

IF (ETA. L£.0.0)ETAa. 000005 
40 ELNaALOG(ETA) 

EKK a 1.386294E0 ♦ ETA * (.9666344E-1 

1 (.3590092E-1 ♦ ETA t (.3742564E-1 

2 .1451196E-1 ))) - ELN * (.5 + ETA 


)/ 2 . 0 ) ) ) 




♦ ETA « 

+ ETA * 

I (.1249859E0 
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onrioonn 


+ ETA t ( .3329355E-1 


f ETA « 


+ ETA « < .5260601E-1 
.173AS06E-1 ) ) ) - ELN 


.52&4496E-2 


+ ETA 
n ) ) 


+ ETA 
t ETA I 


f ««$ S c tC 


3 ETA t ( .6880249E-1 

A .4417870E-2 )))) 

EEK s 1. + ETA 1 (*4432514E0 

1 ( .4757384E-1 + ETA # 

2 <.2499837E0 + ETA * (.9200180E-1 

3 ( .4049498E-1 + ETA 4 

RETURN 
END 

C4»»*»«i**»BEGINNIN6 OF ELEMENT ELLLC/ 

SUBROUTINE ELLC <A#K»E»I) 

THIS SUBROUTINE CALCULATES THE ASSOCIATED COMPLETE ELLIRTIC INTEGRALS 
Or THE FIRST OR SECOND KIND 
THE ARGUMENTS ARE' 

A argument (K SQUARED) FOR UHICH E' OR ' UILL BE FOUND 

K UALUE of ASSOCIATED COMPLETE ELLIPTIC IN'^'EGRAL OF FIRST KIND 

£ VALUE OF ASSOCIATED COMPLETE ELLIPTIC INTEGRAL OF SECOND :.IN 

I IF EQ 1» COMPUTE K i IF EQ 2» COMPUTE E 

DOUBLE PRECISION K i £ > CON ( 32 ) * A » LN4 » CF < 29 ) > CL ( 3 ) » DLOG 
EQUIVALENCE ( CON ( 1) » CF < 1 ) ) » < CON ( 30 ) » CL ( 1 ) ) 

DATA CF /9.AS73590797589018D-2» 3. 088557348 A752694D- 2 » 1.4978983173 
1704629D-2f 9. 65875798 A 1753 11 3D-3» I . 1 2089 13554644092D- 2 • 1 .3855601247 
2l5656D-2»6.6905509906a97936D-3*6 . 499 844 3329390 18D- 4 * 1 .249999999411 
37923D-1 » 7. 03 1242646462736 1D-2M. 88180585654039520-2 #3.7068399934 15 
45422D-2# 2.718986111673825D-2# 1 . 4 105290776153048D-2 » 3 . 1 o3 1 30992 7362 
5896D-3f 1 .50491317836018830-4.4.43147181121559060-1 .5. 6305657874695 
6358D-2»2.1376220647l86198D-2» 1 . 25 1 05924 1 084464 4D-2 » 1 .3034146073731 
7 4 32D-2 » 1 .5377 10252355201 9 0-2 #7. 3356 1649742903650-3. 7.09 9096 409999 7 
322? D-4# 2. 49999999936 17 6220-1 . 9 . 3749920249680 1 1 30-2 » 5 . 6 5S233953655 9 
9024D-2# 4.238280745694790-2. 3.03027477234128430-2 / 

DATA CL /I .5525129948040721D-2# 3. 48336794358964920-3 » 1.642721079 
17048025D-4 / 

LN4 = 1.3862943611198900 
IF (A.EQ.0.0) GO TO 40 
GO TO ( 10 » 20 )>I 

10 K a LN4 + (((((( (CQN(8) *A+C0N(7) >»A+CQN<6> > *A+C0N(5) )tA+CONC 4 >) «A 
.r CON< 3) )»AfC0N<2) >»AfCQN( 1 ) ) *A - DL0G(A>*(0.5P(<(((<(C0N<16)*A+ 
2C0NU5))*A + C0N<14))*A + CQN<13)) *A + CON( 12: ) «A + CON( 11))*A + C0N(10)MA 
if C Q N ( 9 ) ) » A ) 
oQ TO 30 

20 E = I . 000+( < < < ( < <C0N(24 ) *A+C0N<23) )4A + C0.N(22) )<A + C0N(21 ) > *A + CON(20 
I ) )»A+CCN< 19) )*A + CQN< 18) ) »A+CQN< 17) ) »A - DLOG( A)»( (((<<(<C0N(32XA 
2+ CON C 31) ) *A + C0N(30) )*A + CQN(29 ) )*A+C0Nl23) ) «^A4C0N(27) ) 4A + C0N<26) ) * 
3A+C0N(25) )»A) 

30 RETURN 

40 K a 0.99999999D30 
E a 1 .000 
RETURN 
END 

C**$»4t*«$4*Ei£GIN^aNG OF ELEMENT HLAMB/ 

FUNCTION HLAMB (BETA.K) 

C THIS SUBROUTINE CALCULATES THE HEUMAN'S LAMBDA FUNCTION OF BETA AND 
DOUBLE PRECISION A.F.E 
REAL K 

DATA TUOP/0. 6366197724/ 

CALL INEL ( FI #£I #PI . BETA.BETA. 1 .0-K#*2 »0.l#l) 

Aal,0-K »«2 

CALL ELLC (A » F# E# 1 ) 

CALL ELLC <A .F.E.2) 

HLAMB a TUOP*(F*EI +(E-F)«FI) 

RETURN 

END 

C4i4:64*4*S-«*&£GINNING OF ELEMENT INNEL/ 

SUBROUTINE INEL ( F . E # P I . A . PH I , SK I » K3 . K2 . K1 ) 

C THIS SUBROUTINE CALCULATES THE INCOMPLETE ELLIPTIC INTEGRALS OF THE 
C FIRST. SECOND AND THIRD KINDS. THE ARGUMENTS ARE' 


c 

F 

VALUE 

CF 

INCOMPLETE ELLIPTIC 

INTEGRAL 

OF THE FIRST 

KIND 

c 

E 

VALUE 

OF 

incomplete ELLIPTIC 

INTEGRAL 

OF THE SECOND 

KIND 


PI 

VALUE 

OF 

INCOMPLETE ELLIPTIC 

INTEGRAL 

OF THE THIRD 

KIND 

p 

A 

VALUE 

OF 

alpha squared 




r 

PHI 

VALUE 

OF 

PHI 




i; 

SKI 

VALUE 

OF 

K SQUAREDO 




p 

A- ~t 

IF EQ 

0. 

DO NOT COMPUTE PI i 

IF NE 0. 

COMPUTE PI 


■j 

K2 

IF EQ 

Of 

DO NOT COMPUTE E i 

IF NE 0. 

CO.MPUTE E 


c 

M 

IF EQ 

0# 

no NOT COMPUTE F ; 

IF NE 0. 

COMPUTE F 



DOUBLE 

PRECISION 

ARG.FDfED 





DATA PIT/1.57079633/ 

£ = 0.0 

R=0.0 

IF (K3.EQ.0) GO TO 10 
CALL ELINT3 < SK I # - A » PH I . P I ) 
iO IF (M .EQ.O) GO TO 30 

IF < AB3<PHI-PI T) .GT. 10 .0»*( -7 ) > GO TO 20 
•»RG=-1 .0-SKI 

CALL ELLC ( ARG . F D . ED » 1 ) 

FaFD 

30 TO 30 

20 CALL ELIHT3 < SK I » 0 . 0 . PH I # F ) 
iO IF (K2.EQ.0) GO TO 50 

Ip (ABS< PHI-PIT ) .GT . 10.0*« ( -7> ) GO TO 4Q 
ARG=1 .0-SKI 

CALL ELLC ( ARG » FD » ED # 2 ) 

C=ED 

■JO TO 50 

40 CALL ELINT3 ( SK I . - 3K I . PH I . E ) 

Ea(1.0-SKI)«E + 0.5»SM*SIN<2.0*PHI) /SORT* 1 . 0-SK I *S I N < PM I > * ♦ 2 ) 
30 RETURN 
END 


100 




:«<i»<<ftt<x&EGlNNlNG OF ELEHENT MATX/ 

SUBROUTINE MATRIX < XPKT f YPKT > 

5 * COMPUTE MATRIX AtB.Z OR X*Y.Z 

COMMON /BLOCKl/ 1 GEOMF < 9 ) . I S I GF ( 9 ) 

COMMON /BLQCK3/NQAXI f NQCftOS . NQUQR T » NOV I * NQV2 * F IftSTE * CASTE * LSSE * SO # 
1 NSMALL 


COMMON/FLG/ 
LrFLGOS 
[ f FLGOB ' 

;»FLG13 / 

. • FLG 1 8 


HEDRM 15) 
FLG04 
FLG09 
FLG14 
FLG19 


»FCG23 »FLG24 

COMMON/NNN/ NTi 
INERlf NERO I 

INUNC(5)» TYPEC<5)» 
1TYPEEC(5) »NUNEC(5) 
DOUBLE PRECISION HEDR 


( CASE 
.FLG05 
»FLG10 
•FLG15 
. FCG20 
*FLG25 
Nil ( 1 1 ) » 
NMA» 

NCF(ll) f 


t NB 

» FLQ06 
iFLGll 
» FL G 1 5 
•FLG21 
»FLG26 
MN » 

NSIGAt 

lECi 


»NNU 
rFCGO? 
» FLOl 2 
*FCG17 
* FLG22 
»FLG27 
NUNA(5 ) , 
NSIGC* 
NSIGCCf 


TYPEA(5) » 


INTEGER 

FLG03 

♦ FLG04 

»FLG05 

.FLGO& 

.FLG07 

MTRX 

240 

1 

f FLG08 

f FLQ09 

»FLG10 

.FLGll 

fFLGl2 

MTRX 

250 


»FLG13 

»FLG14 

• FLG15 

.FLGIA 

»FLG17 

MTRX 

260 

3 

»FLG18 

.FLG19 

»FLG20 

.FLG21 

»FL622 

MTRX 

270 


»FLG23 

»FLG24 

»FLG25 

.FLG26 

»FLG27 


REAL MN 

logical NOAXI iNOCRQS! 
REAL LSSE 


NQVORT* NQUl * NQU2 » F IRSTE * CASTE t PF 


COMMON /TL/ 


COMMON 


A(500) » 3(500) * 

CX<500) » CY(500) 

If J » 

DXf DYf 

XK » ££K » 

/KIT2EL/ UXA(5> »UYA<5) » 

yXC(5) ,VYC<5) |02C(5) 
VN<5) »UT<5) 


AX (500) » 
CZ(500> » 

J1 » 

NI . 

£KK f 

S1GA(500f5) » 
»SIGC(S00»5) 


Ar(500) 
AXV( 

SJ» 

XJ. 

K 


AZ(500) » 
00) » AYU(500> f 
DS» 

Y J* 

PF 


« START 
* INITIALIZE 
L.O 120 J » 1 »5 
VN(J) » 0»0 
120 UT(J) * 0. 

J1*0 
N1 » 0 

initialize the entire rou 
DO ISO J»1»NT 
AX(J) a 0.0 
AY( J) a 0.0 
A2(J> a 0.0 
CX<J) a 0.0 
CY( J) a 0.0 


CZ< J) 

a 0. 

0 

AXV( J) 

a 0 

.0 

Arv< J) 

a 0 

.0 

DO 220 

K=l 

>NB 


1) N0U2 a .TRUE. 

2> NuOl a .TRUE. 

THE columns of THIS ITH RQU aS INFLUENCED 


- 1 


C» 


NOUl a .false. 

N0V2 a .FALSE. 

IF (ISICF(K) .EO 
IF (ISIoF(K) ‘cQ 
BEGIN MATRIX FORMULATION OF 
BY THIS MH BODY .. 

Mi a N1 F 1 
N1 = N1 + ND(.<) 
rIRSTE .TRUE. 
lASTE a .FALSE. 

DO 170 JaMlfNl 
= Ji + 1 

IF (J .£Q. Nl) LASTE 
CALL PAKAB(XPKT»YPKT) 

170 FIRSTE a .false. 

J1 a Ui + 1 

MATRIX FORMULATION OF COLUMNS FOR THIS BODY HAVE BEEN COMPLETED. 
IF(FLGl5.LE.0>G0 TO 220 

:f (nlf(K) .gt.o) go to 220 

DO 190 J a Ml» Nl 

UN(K> a UN(K)+AXV(J) 

190 VT(K) a VT<K)fAYV(J) 

220 CONTINUE 
RETURN 
END 


TRUE . 


*«$»<<tC#BEGINNING OF ELEMENT PRAB/ ttttttttttC 

SUBROUTINE PARAS ( XPKT » rPM ) 

OVERLAJ' (PAN»1»2) 

PROGRAM PARA&(XPKTfYPKT ) 

CONTROL FOR PARABOLIC X»YfZ MATRICES COMPUTATION. 

COMMON 7CL/ Xl(500)» Yi<500)» X2(500)r Y2(500>. BELL<500). 

1 SINA(500> .C0SA(500) 

COMMON ,/TL/ A(500)» B<500>» AX(500)» 

1 CX(SOO)» CY(500)f CZ(500)» 

3 If Jf Ji. 

4 D X ♦ [I Y ♦ N I » 

XK f EEK » EKK » 


AY(500)f AZ(SOO)* 
AXV(500> fAYVCSOO) » 
SJ* OSf 

XJ» YJ» 

K » PF 


COMMON /BL0CK2/ HCURV ( 500 )» HARD ( 500 ) 
common /BLQCK3/N0AX I t NQCROS . NOUOR T f NGVl » N0U2 » F IRSTE » CASTE » LSSE » 50 » 
INSMALL 

logical NOAXI»NOCROS»NOVORT»NOU1 f N0U2 » FIRSTE fLASTE 
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REAL 
C OFF-BODY 
JlPl 


LSSE 


Jl + I 

D1 - ( XPKT-Xl < J1 ) ) »*2 + (YPKT-Y1(J1>)41*2 
D2 s (XPKT-X2( J> )»*2 + ( YPKT-Y2 < J ) ) ««2 
D3 = (XPKT-XllJlPl) )»*2 +• (YPKT-YI(J1P1>>**2 
PMINSQ * D1 

IF (D2 .LT. DrtINSQ) DrtINSO - D2 
IF (D3 .LT. DMINSO) DttINSQ = D3 
OMIN =s SQRKDHINSQ) 

•F (DhlN .EQ. 0.0) GO TO 60 
NI = 16.*HARC( J)/DrtIN + 0.9 
IF <NI .GT. 0) GO TO 50 
NI - 3 
D3 = HARC(J) 

GO TO 80 
50 NI s 2*NI 


IF 
60 NI 


90 XJ * 
YJ = 


(NI .LT. 128) GO TO 70 
= 129 

OS = HARCCJ)/64. 

GO TO 80 
70 XNI = NI 

DS * 2.«HARC< J)/XNI 
NI - NI + 1 
X2< J) 

Y2< J) 

30 = -HARC(J) 

CALL PARAB2(XPKTiYPKT) 

RETURN 

END 

C**»S|!»«1»»BEGINNING of element PRAB2/ 

SUBROUTINE PARAP2 ( XPK T , YPKT ) 

C THIS ROUTINE CALCULATES THE UCLOCITY CONTRIBUTIONS OF 

C OFF-DIAGONAL (AND ENDS OF THE ON-DIAGONAL) ELEMENTS BY A 

C SIrtPSON-RULE. ELEMENTS MAY BE FLAT OR CURV£D» WITH PIECEWISE 

C CONSTANT-* LINEAR-* OR QU ADR A T IC -SOURCE DENSITY. 

COMMON /BLOCK2/HCURU(500 ) » HARC(SOO) 

COMMON /BL0CK3/NQAXI * NQCROS » NQVOR T , NOV 1 * N0V2 » F IKSTE * LASTE » LSSE » 30 » 
1 NSMALL 

..OMMON /BLOCK4/A11(500)*A12(500) * A13(500) » A21 (500) * A22(500) » 

• A23(500) 


COMMON /CL/ 

XI (500) * 
3INA<500) 

Y1 (500) , 
*CQSA<500) 

X2(500) * 

Y2(500)» D£LL(500 

‘common /TL/ 

A(500) * 

B ( 500 ) * 

AX(SOO) » 

AY(500)* A2(500>* 


CX(500) * 

CY(500)* 

C2(500) * 

AXV(500) *AYV(500) . 

2 

I * 

J * 

Jl * 

SJ* DS* 


DX* 

DY* 

NI* 

XJ* YJ* 

J 

XK 1 

EEK * 

EKK » 

K* PF 


LOGICAL HOAX I * NOCftOS * NOVORT * NQVl * N0U2 » F IRSTE » LASTS 

logical PF, ysmall 
PEAL LSSE 

DIMENSION OOAXIX(129)*UOAXIY(129) » VOCftSX ( 129 ) *OOCRSY( 129) * 

1 U0CRS2 ( 129 ) * U0O0RX< 1 2 9 ) * ‘.'OVOft Y < 1 29 ) * XR I NO < 1 29 ) » YR INC ( 1 29 ) 
DATA PI/3.1415927/ 

OFF-BODY ... 

X = XPKT 
Y a YPKT 


5INALaSINA( J) 

CQSAL=COSA( J) 

HC = HCURO(J) 

'if (FIRSTE) Msj+l 
IF (LASTE) rt=J-i 
mMl*M-l 
MP1=M+1 

C BEGIN BY GENERATING THE INDUCED VELOCITIES AT THE ENDS 3^ Ml 

C INTERVALS (OF UNIFORM LENGTH) 

3 * SO 

DO 110 L = 1»NI 

ssa=s *#2 

C CALCULATE LOCAL SOURCE RADIUS* QA* AND AXIAL LOCATION* OB. 

QA = YJ + SINAL»S+HC*COSAL*SSQ 

jB = XJ + COSAL*S-HC*SINAL«SSa 

XRING(L) = OB 

;'RIN6(L) = QA 

IF (OA .GT. 0.0) GO TO 39 

vOAXIX(L) - 0.0 

VOAXIY(L) = 0.0 

VOCRSX(L) = 0.0 

VOCRSY(L) * 0.0 

VOCRSZ(L) == 0.0 

VOVQRX(L) = 0.0 

VOVORY(L) = 0.0 

GO TO 110 

C CALCULATE THE COMPLETE ELLIPTIC INTEGRALS* K(K> AND E ( K ) * 

C OF THE FIRST AND SECOND MNDS.... 

39 XMB=<-QB 

XMBSQ= XMB*«2 

YMASQ = (Y-QA)**2 

rPASQ*(Y4-QA>*»2 

n 3 YPASQ + XMBSQ 

XK= 4.*0A«Y/T1 

■CALL £LIP(XK*EEK*EKK) 

C SUBROUTINE SLIP RETURNS THE VALUES AS EKK AND EEK . 

rSQ = Y*»2 
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ASQ 3 QA<«2 
ROOTl » SQRT(Tl) 

T2 a rMASQ > XnBSQ 
T3 = XMBSQ + ASQ 
ai (-4.4tEEK)/<R00Tl*T2) 

Q2 = <2./R00TU*(EKKf<YSQ-ASQ-XHBSQ>*€EK/T2) 
rSMALL = CrSQ/T3) .LT. .0001 
IF (YSMALL) GO TO 40 
■JO TO 50 

ROOT32 = (S0RT<T3) )»»3 
T4 = FI«QA/RQ0T32 
T3SQ s T3**2 

IF (NOAXI) go TO 70 
IF (YartALL) GO TO 60 
OOAXIX(L) s Ql»OA*XrtP 


Ql»OA*XrtP 
■v'OAXlY(L) - -Q2»QA/Y 
JO TO 70 

OOAXIXIL) - l-2.*XhP*T4)*<1.0 + 0. 75» r SO* < 3 . 0 #hSQ- 2 . *XH&S3 )/ T35G > 
OOAXIY(L) = -Y*T4*(2.fXMBSQ--ASQ)/T3 
IF (NOCRQS) GO TO 90 
IF (YSMALL) GO TO 80 

OOCRSX(L) = 2.*(XMB / ( Y»RQQT 1 ) ) * i EKK- ( YSQ+T3 > *EEK/ T2 ) 

V0CRSY(L>*2. *< 1 ./< YSQtROQTl > ) *< -( T3>*£KKi( T J*»2f YSQ* ( XM&SQ-ASQ > ) 
1 4EEK/T2) 

00CRSZ(L)»2.* (ROQTl/YSa)«< ( YSQ T3 ) *£KK/T 1 -EEK ) 

00 TO 90 

OOCRSX(L) * -3 . f QA*XMP*Y1!T4/T3 

'JOCRSY(L) = QA*T4*ll.0-0. I 25 » Y SQ * t 4 . *XMB3Q -ASQ > / T3SQ ) 

=JOCRSZCL) s QA«T4«(1 , 0 + 0 . 375 4 YSQ4 (. ASQ-4 . *XhPSQ ) / T3SQ ) 

IF (NOOQRT) GQ TO 110 
'^OOORX(L> * (ASQ-YSQ)*Q1 - 02 
IF (YSMALL) GQ TO 100 
OOOORY(L) s Y»XMB*Q1 + XMB*Q2/Y 
GQ TQ 110 
OOyORY(L) a 0.0 
S = S + DS 

assemble THE ROW OF EACH DESIRED ttATRIX... 

IF (NOAXI) GQ TO 120 

CALL SIMS0N<V0AXIX»NlfDSfU0AX) 

CALL S!MSQN<OOAXIYiNIfDS»OOAY ) 

AX(J) = AX<J) + VOAX 

AY(J) = AY(J) + OOAY 

IF (NOCROS ) GO TO 130 

CALL SIMSQN(VOCRSX»Nl»DSiOOCX) 

CALL SIMSON<OOCRSY»NI»DS»VOCY) 

CALL SIrtSON< V0CRSZ»NI »DS» V0C2) 

CX(J) » CX<J) + VOCX 

CY< J) = CY< J) + OOCY 

CZ( J) » CZ( J) + 00C2 

IF (NOOORT) GO TQ 14U 

CALL ■ S1MS0N( VOVORX »NI > DS * OOVX ) 

. CALL SIMSON( VOOORYfNI »DSrOOOY) 

AXO(J) = AXO<J) + OOVX 
HfO(J) = AYO<J) + OOvl 
CONTINUE 

IF (NOUl) GO TO 250 
IF (NOAXI) GO TO 170 
S - SO 

DO 160 L a liNI 
OOAXIX(L) = VOAXIX(L)XS 
/OAXIY(L) a V0AXIY(L)«S 
S a S + OS 

CALL SIMSON ( UQAX I X » N I ♦ DS » 01 AX ) 

CALL SIMSON ( VOAXIY .NI » DS » UlAY ) 

AX(MMl) a AX<Mttl) + A11(J)*U1AX 

AY(MMl) a AY(MMl) + All(J)*OlAY 

•^X(M) a aX(M) + A12(J)«U1AX 

AY(M) a AY(M) + A12(J)»01AY 

AX(MPl) a AX(MPI) + A13(J)*U1AX 

AY(MPi) a AY(MPI) + A13(JMU1AY 

IF (NOCROS) GQ TO 200 
a SO 

HO 180 L a 1 ,NI 
vOCRSX(L) a V0CRSX(L)»S 
.»OCRSY<L) a OOCRSY(L>»S 
vOCR3Z(L; a OOCR3Z(L)»S 
S a s + DS 

CALL 31MS0N(V0CRSXfNI »DS#V1CX> 

CALL SIMSON(OOCRSY»NIfDSfUlCY) 

CALL £IMS0N( VOCRSZ»NI »D3» VICZ) 

CX(MMl) a CX(MMl) + All(J)»UtCX 

CY(MMl) a CY(MMl) + A11<J)«01CY 

CZ'MMl) a CZ<MM1) + AIKJMVICZ 

CX(M; a CX(M) + A12(J)»01CX 

CY'M) a CY<M) + A12(JXU1CY 

CZ(M) a CZ(M) + A12(J)*01CZ 

CX(MPl) a CX(MPl) + A13(J)*V1CX 

Cr(rtPl) a CY(MPl) + A13(J)4V1CY 

■:Z(MP1) a CZ(MPl) + Ali<J)«VlC2 

IF (N0U2) GO TO 2S0 

IF (NOAXI) GO TO 220 

S a so 

no 210 L a i,NI 
VOAXIX(L) a V0AXIX(L)*S 
OOAXIY(L) a UOAXIY(L) fSOO 
3 a 3 + D 3 

CALL SIMSON ( UOAX I X . N I » DS • 02 AX ) 

CALL SIMSON ( VOAX I Y p N I » DS » U2AT ) 
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AX(MMl) • AX(MMl) + A21(J)#»J2AX 
AY<MMl) 3 AY<HHl) t A21(J)*W^2AY 
AX(H> 3 aX(H) t A22<J)tg2AX 
AY(rt) 3 aY(«) + A22(J)»U2AY 
AX(MPl) 3 AX(MPl) + A23<J)«U2AX 
AY(rtPl) 3 AY(rtPl)f A23(J)<y2AY 
220 IF (NOCROS) 60 TO 250 
3 3 SO 

DO 230 L 3 1 ,NI 
VOCRSX<L) = yoCRSX(L)*S 
OOCRSY(L) 3 OOCRSY(L)*S 
00CR3Z(L> 3 UOCRSZ(L)*S 
230 S 3 S + DS 

CALL SIMSON< UOCRSX »NI » DSiOZCX) 
CALL SlMS0N(V0CRSY.NlfDSig2CY) 
CALL SlhS0N(00CRSZiNIiDSt02CZ) 
CX<Hhl) 3 CX(rthl) + A21(J)*V2CX 

CY<rtrtl) 3 CY(ttttl) V A21(J)»V2CY 

CZ(MHl) 3 CZ(MMl) + A21(J)*g2CZ 

CX(rt) 3 CX<M) + A22<J)«02CX 

CY(M) 3 CY(M) + A22(J)*g2CY 

CZ(M) 3 CZ(rt) + A22(J)*V2CZ 

CX(MPl) * CX(MPl) + A23<J)*g2CX 
CY(MPl) 3 CY(HPl) 


C2(MP1 ) 3 C2(MP1 ) 


A23( J)*y2CY 
A23< J)*V2CZ 


250 


RETURN 
END 

C»*$$»*$$Y*&EGINNING of ELEMENT QCC/ 
SUBROUTINE QC ( QMEG f QM i Q ) 




C THIS SUBROUTINE CALCULATES THE LEGENDRE FUNCTIONS OF THE SECOND KIND 
C AND HALF ORDER. THE ARGUMENTS ARE' 

C OMEG ARGUMENT FOR WHICH LEGENDRE FUNCTIONS WILL BE FOUND 

C QM VALUE OF LEGENDRE FUNCTION OF MINUS ONE HALF ORDER 

C Q VALUE OF LEGENDRE FUNCTION OF PLUS ONE HALF ORDER 

DOUBLE PRECISION OMEGD « ARG f A * F « £ « QHD > QD 
GMEGD=OM£G 
ARG = 2.0/(0MEGD-K.0) 

A31 ,0-ARG 

CALL ELLC <A»F»£»l) 

CALL ELLC (A»F,E»2) 

•nMD = F»ARG»*0.5 

JD=-E*< 2. 0*<0M£GD+1 .0) ) » yO . 5+QMEGD*QMD 

QM3QHD 

a = QD 

RETURN 

END 

C#*$$*««$»*&EGINNINQ OF ELEMENT SIMSQ/ 

SUBROUTINE S IMSON ( Y » N » DX . ARE A ) 

C THIS ROUTINE INTEGRATES Y OVER \N-l) INTERVALS OF EQUAL LENGTH» DX» 

C YIELDING THE ENCLOSED AREA. (N MUST BE An ODD INTEGER.) 


DIMENSION Y(l) 

C 

SUM 3 Y<i) + 4.0«Y(N-1) + Y<N) 

IF (N .EQ. 3) GO TO 20 
NM2 3 H - 2 
DO 10 I=2#NM2»2 

10 SUM 3 SUM + 4.0*Y<I) + 2.0*Y(I-H) 

C 

20 AREA 3 < aBS(DX)/3.0) * SUM 
RETURN 
END 

SUBROUTINE LINEAR< X» Y » IMAX» IMIN) 

DIMENSION X< 1) » Y< 1) 

IM=IMAX-1 
NTOT=IMAX-IMIN 
URITE(6»199) NTOT 
199 FORMAK 15) 

DO 10 I3IMIN»IM 

XMN = 0 .5»<X< I + l ) -X< I) ) 

YMN30. 5*( Y< I + l )-Y< I ) ) 

AX = ABS< XMN) 

IF(AX. LE.O. 0000001 ) GO TO 10 

rp:-YMN/XMN 

XMN*X< I ) +XMN 

YMN=Y(I)+YMN 

WRIT£(6»93) I » YMNf XMNf YP 
93 F0RMAT<2Xf I5#3(lPEi2.4) ) 

10 CONTINUE 
RETURN 
END 

SUBROUTINE EFFICY ( NCOEF r NP TS » NS ) 

IMPLICIT R£AL*8(A-HfO-Z) 

REAL# 4 XINf YYP»YYD»YOUT»YPOUT»XrtAX»xr:N»x(lAfXMIf DF D»PLWC»£PT» 
1 XA< 1000) # YA< 1000 ) 

REALX4 XBTA(iOO) »YBTA<100) 

COMMON/COEF/A(20»21 ) 

COMMON/COEF1/B120f 30) 

C0MM0N/C0EF2/XMA( 10) »XMI < 10) 

COMrtON/DPDl/NSIfPLUC( 1 0 ) f I R3 ( 1 0 > f DPD ( 1 0 ) r EP T 
REWIND 3 


Nl3i 

N2 = 0 

NN3NC0EF+1 
DO 1 IJK«l»NSI 
NN13NC0EF+IJK 

IFdJK .£Q. NSl/2) NDA=IRS(IJK) 
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ft£AD(3) NOATA 

IF(NS.e0.1 .OR. NSr .GT. 1) NPTSaNDATA 
N2»NDATA-fN2 
00 543 I«N1»N2 
R£A0(3) XAU ) f YA( I) 

rF(APS(XA(D) ,L£ .0.0000001 ) XA( I > »0.0 
543 CONTINUE 

IF<NSI.GT.n URITE<6»72)IJK 

IF(NSI .GT. 1) yRITE<6r 71) DPD( UK) »PLUC( UK) 

URITE(6»73) NDATA 
UR1TE(6»90) 

00 97 I=N1.N2 

97 URITE<6»93) I » YA < I ) f X A ( I ) 

IF(K‘SI .EQ. 1) 60 TO 2 
XhIN=XA(Nl > 

XrtI < UK ) sXA(Nl ) 

XHAX*XA(N2) 

XMA<IJK)*XA<N2) 

IF(XA(Nl>.LE.XrtIN> XhIN=XA<Nl) 

IF(XA(N2) .GE.XMAX) XMAX=XACN2) 

XIN = XA< N1 ) 

CALL TERP(XINf YOUT» YPOUTf XA» YAf N1 »N2»NC0EF»NPTS) 

DO 3 II=1»NCQEF 
B( II iNNl )=A< II»NN) 

3 CONTINUE 
N1«NDATA+N1 

1 CONTINUE 

URITE(6»198) XMINfXMAXiEPT 

CALL EF<XHAX»XHIN.NDA*NSI#NCQ£F»PLWC) 

RETURN 

2 URIT£<A»99) 

URITE< 6i 91) 

NT0T=N2-N1+1 

CALL LINEAR ( XA > YA » N2 i N 1 ) 

URIT£(A*70) NCOEFjNPTS 

IAC = 0 

NTQA*NTOT 

NN-NCOEF+1 

URITE<4#91) 

DO 20 I=NX»N2 
IAC=IAC+1 
XIN=XA< I ) 

IF( NS.EQ . 1 .AND. I .GT. 1) GO TO 21 

CALL TERP( XIN» YQUTi YPOUT »XA# YAf N1 f N2 r NCOEF r NPTS ) 

21 IF(NS'.£0.1) GO TO 22 
YYP-YQUT 
YYD*-YPOUT 

GO TO 23 

22 YYD*0.0 
YYPaA< 1 »NN ) 

DO 9 Js2»NC0EF 

IF<XA<I) .EO. 0.0) XA< I)sO. 0000001 
YYP»YYP+A(J#NN)4XA<I)4»<J-1) 

YYD*YYD+ ( J-l )#A<J»NN)*XA<I)*f<J-2) 

•? CONTINUE 
YYD*-1 . tYYD 

23 IF<YTD.LT .0.0) GO TO 27 
yRITE<6.93) I f YYP» XA( I) f YYO 
GO TO 28 

27 NT0A=NT0A-1 
tAC=IAC-l 

IF( I .EQ . N1 . OR . I .cO. N2) GO TO 29 
URITE<A»93) I » YYP»XA( I ) f YYD 
GO TO 20 
29 YYD^O.O 

yRITE(6»93) I » YYP»XA( I ) > YYD 

NT0A=NT0A+1 

IAC=IAC+1 

28 YBTAdAOaYYD 
XBTAC IAC)=XA( I ) 

20 CONTINUE 
URITE<6>92) 

URITE(6f 199)NT0A 
URITE<9)NT0A 
DO 24 IslrNTOA 

URITE<6»95) I »XBTA< I ) »YBTA< I ) 

■JRITE<9) XBTA< I) » YBTA( I ) 

24 CONTINUE 

C FORhATS 

70 F0RNAT(/»2Xi 'POLYNONINAL LEAST SQUARE FIT'f//*2X* 

1 'NO. OF COEFFICIENCTS IS 'tlZf/f 

2 2X»'N0. OF POINTS FOR DATA FITTING IS'fI3»/> 

71 FORMAT </r2Xf 'DIAMETER' f 4 X» 'PERCENTAGE' f 

1 /»2Xf 'DISTRIBUT. ' f2Xf 'OF U. CONT.'f 

2 //f2Xf2<1PE10.3f3X) ^/) 

FORMAT(///»lXf '»*»FOLLOUING OUTPUT IS FOR DROPLET SIZE 

DISTRIBUTION — ' f I3f/ » 

3X»'(THE DATA ARE USED TO DO CURUE FITTING)'./) 
F0RMAT(/f2Xf 'NO. OF DATA POINTS STORED IN lAPEOS » 'fI6 
F0RMATC/f5X» 'I ' f6Xf ' ZO '»7Xf' SD'f/) 

FORMAK / f 5X F ' I ' r6Xf ' ZO 'f7Xf' SD'fSXf' BETA'f/) 

FORMAT (/f5X f ' I ' f 6X f ' SD ' f 8X f ' BETA ' f / ) 

FORMA T ( 2Xf I5f3< 1PE12. 4) ) 

F0RMAT<2X.I5»2(1PE12.4) ) 

F0RMAT</i2Xf 'LINEAR APPROACH'f/) 

F0RMAT(///»1Xf 'FOLLOWING OUTPUT IS FOR THE O'/ERALU 
1 MULTI-SIZE CASE' »//f4Xf 'UPPER SURFACE LIMIT'f?Xf 
2'LOWEft SURFACE L I MI T ' f / f 1 2X f ' SU ' . 2AX f ' SL ' f / f AX f 
32< 1PE12.4 . 17X) f//. 2X# ' TOTAL COLLECTION EFFICIENCY'f 


72 


73 

90 

91 

92 

93 
95 
99 

198 
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43X*lPE12.S*/> 

199 FORMAT< 15) 

RETURN 

END 

SUBROUTINE EF ( XMAX f XM I N » NOA . NS I »NCOEF»PLUC) 

IMPLICIT REAL»8( A-HfO-Z) 

DIMENSION XA( 100) »PL'JC( 1 ) . XPTA( 100) ♦ YPTA< 100) 

REAL«4 XMAXf XMIN#PLUCf YYDtDEStXMAtXMl »XA 
COMMQN/CO£F1/B(20» 30 ) 

C0MMQN/C0EF2/XMA( 10 ) • XMI < 10 ) 

DES*( XMAX-XMIN) /FLOAT (NBA- 1 ) 

IAC = 0 

ntoa«nda 

URITE(6»70)NCCEF 

URITE(6»92) 

DO I IJ*1»NDA 
lACalAC+1 
rYD = 0.0 

XA( IJ)*XMAX-PES*(NDA-IJ) 

DO 2 I=1»NSI 

IFCXACIJ) .GT. XMACD.OR. XA(lJ).Lr. XMKD) GO TO 2 

nn=ncoef+i 

DO 2 J*2»NC0EF 

1F(XA(IJ) »EQ. 0.0) XA(IJ)sO. 0000001 
YYD=YYD+( J-1)*B( JfNN)*XA(IJ)»*< J~2)*PLUC<I) 

2 continue 

YYD=-1 .*YYD 

IF<YYD .LT. 0.0) GO TO 27 
URITE(A»93) lJfXA(IJ)fYYD 
GO TO 28 
27 NT0A=NT0A-1 
IAC=IAC-1 

IF<IJ.EQ.1 .QR.IJ.EO.NDA) GO TO 29 
URITE(4»93)IJ#XA(IJ)»YYD 
GO TO 1 
29 YYD=0.0 

URITE(6»93) IJ»XA(IJ)»YYD 
NTQAaNTOA+l 
lAC*IAC+i 
29 YBTA(IAC)=YYD 

XBTA( IAC)=XA( IJ) 

• CONTINUE 
'ORITEt 6»92) 

URITE< Af 199)NT0A 
URITE(9)NT0A 


24 

:o 


^2 

?3 

199 


C 


DO 24 Isl»NTQA 

URITE(6»95)I«XBTA<I)»YBTA<I) 

URITE(9) XBTAa ) » YBTA< I ) 

CONTINUE 

FORMATS 

F0RMAT(/»2X» 'POLYNOMIAL LEAST SQUARE FIT'*//»2X* 

'NO. OF COEFF • CIENCTS IS '»I3»/» 

\ 2X»'UITH WHOLE DATA SETS'»/> 

FORMAT</ f 5X » ' I ' »AX » ' 3D ' f 8X » ' BETA ' »/ ) 
format (2X» I5»3( 1PE12.4) ) 

FORMAT <2XfI5»2<lPE12.4>> 
format ( 15) 

RETURN 

END 

SUBROUTINE TERP ( X IN » YOUT » YPOUT , X A , Y A f N I r N2 » NCOEF t NPTS > 

DOES simple NCOEF POLYNOMIAL LEAST SQUARES FIT TO NPT3 OF XA» 
IMPLICIT REAL*8<A-H,0-Z) 

REAL *4 XINf YOUr f YPOUr^XA( 1 ) r YA( 1 ) 

DIMENSION F( 1000»20) » Y ( 1000) 

C0MMON/C0EF/A<20»2l ) 

NCMI s NCOEF - 1 
NCOL s NCOEF f 1 
SEARCH XA FOR XIN 
DO 1 IF»NlfN2 


I^IF 


iF(XIN.LE.XAdF) )G0 TO 5 

CONTINUE 

CONTINUE 

IMIN«I-NPTS/2 

iFdMIN .LT. Nl) IMIN^Nl 

IMAX«IMIN+NPTS-1 

IF< IMAX .GT ,N2) IMAX=N2 

IMIN«IMAX-NPTS+1 

III « 0 

DO 211= IMINp IMAX 
III = III + 1 
FdII»l) = 1.0 
DO 3 JJ = 2fNC0EF 
FdllpJJ) = XAdI)»*(JJ-l) 
Y d II ) * YAd I ) 

A( 1 p I ) = NPTS 
DO 10 J * IpNCOEF 


DO 20 I » JpNCOEF 

iFd .EQ. 1 . AND. J.EQ. 1 )G0 TO 20 

SUM « 0.0 

DO 11 ISUH * IpNPTS 
:i SUM - SUM + FdSUMpIXFdSUMpJ) 
A( Jp I ) = SUM 
Ad p J) = SUM 
20 CONTINUE 
SUM ® 0.0 
DO 50 I = 1 pNPTS 
^0 SUM = SUM + Y<I)#FdtJ) 
A<J»NCOL) » SUM 


YA 
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iw CuNFIML'E. 

CALL CHOLt5(A»rtCOEF»NCOL»l) 

YOUT*A(IiNCOL> 

YP0UT*0.0 
DO 4 I=2»NC0EF 

IF<XIN .£0.0.0) XINsO.OOOOOOV 
YOUTsYOUT+A(I »NC 0L) *XIN**( I-l ) 

4 YPOUT»YPOUT + ( I-l)*A<IiNC0L)*xnu*<I-2) 

return 

END 

SUBROUTINE CHOLES( AtNrHfhATSYn) 

CDPD MATRIX SOLUTION ROUTINE FROM DON TO£»f» ON 4-24-00 

REALMS A(20»21) »SUMi temp 
IF < A( I , I ) .NE .0 .0) GO TO AV 
DO 37 J=2»N 

IF<A( J,l) .EQ.O.O) GO TO 37 
IFLIP=J 
60 TO 27 
37 CONTINUE 
GO TO 54321 
27 DO 57 K=l»tt 

TEnP=A( IFLIP»K) 

A( IFLIP »K )*A< 1 »K) 

A( I )=TEMP 
57 CONTINUE 
47 DO 2 J=2»rt 

A( 1 > J)-A( 1 1 J)/A(l I 1 ) 

2 CONTINUE 
DO 6 Ia2»N 
DO 7 Ja2fH 

IF(MATSYf1,EQ.0)G0 TO 49 
IF< I-J)69»4a»67 
49 IF< J.CT.DGQ TO 69 

68 K=J-1 
SUM^-0 ,0 

DO 3 IftsliK 

SUn-SUhiA( I f I R) «A ( IR f J ) 

3 CONTINUE 

A<lt j)aA(I» JJ-SUH 
GO To 7 

69 K=I-1 
SUrt*0 . 0 

DO 4 IRsliK 

SUh=SUn-VA< IilR)*A<Ift»J> 

4 CONTINUE 

IF ( A( I , I ) .EQ . 0 . 0 ) GO TO 54321 
A(I»J)=<A<IfJ)-SUM)/A(IiI) 

GO TO 7 

67 A ( I » J) sA( J • I ) «A ( J I J) 

7 continue 

6 continue 

DO 52 K»2»N 
I=N+1-K 
SUri = 0.0 
LL*I+l 

DO 51 IR=LL»N 

£Ui1sSUn + A( IiIR)*A<IR»M) 

51 continue 

A(I »H> -A( 1 irt ) -SUfI 

52 continue 
GO TO 12345 
N»-l 

return 

END 

SUBROUriNE APLOT ( U # YL » N » h » LDR » FC # C2 ) 

DIMENSION U ( I ) r YL < 1) »N< 1) 
dimension xri TL< 20) » YTITL(20) 

DIMENSION DSNAME<9) 

DATA IDP/1/»DSNAME/'SAUE' » 'PLOT' »7»' '/ 

K: no. of DATA SETS 

N(K): NO. OF DATA POINTS IN KTH DATA SET 
U: X-COORDIHATE 

yl: y-coordinate 

1F(LDP .EQ. 1) GO TO 19 
READ(2flOO) UTITLa)»I*l»20) 

READ<2fl00) ( YTITL< I ) » I=l » 20) 

REAlM 2f 101 )YHEIT » YhAX i YMIN 
READ(2 I 101 ) XLENGf XhAX f XMIN 

100 FURMAT(20A4) 

101 FORMAT ( 5< 5X fFlO .0 ) ) 

C2a(YhAX-YMlN)/YHEIT 
FC* t XMAX-XMIN )/XLENG 
CALL CCD£GN( lO.Of 1 .0) 

CALL BEGIDCIDP)^ 

CALL AXIS(OiIoiVxTITLf-80iXLENG»0.»XMIN.FCfO.O) 

CALL PL0T(0. »0. »3) 

CALL AXIS<0.»0.»YTITL»80»YHEIT»90.»YMIN»CZ»0.0) 

CALL PL0T(0. »YHEITf3) 

CALL PL0T(XLENG»YHEITi2) 

CALL PLOTlXLENGfO. t 2) 

CALL PLOT(O.0»0.Of3) 

19 continue 
NK s- I 
JN-0 

DO 14 I «1 » K 
Jl=JNf 1 
JN« JN+N( I ) 

NK»NK4-2 

DO 13 JsJliJN 

2=( Yl(J)-YMIN)/CZ 

X=(U( J)-XMIN)/FC 

IF(J .EQ. Jl) call PL0T(XiZ»3) 

IF(Z.GT.YHEIT) GO TO 13 
IFlLDP.EQ.l) 00 TO 17 
IF ( I .NE . K) GO TO 67 
17 CALL PL0T(XiZ»2) 

UO TO 13 
67 continue 

call SYMB0L(XiZ»0.14»NKf0. »-l ) 

13 continue 

CALL PL0T(0.0»0.0»3) 

14 CONTINUE 

CALL ENDID(IDPi“1»DSNAME) 

CALL CCBECN( -1 .0 » 1 .0) >1 

RETURN IU7 

END 

SUBROUTINE TRAJCT 

return 

END 


34321 

12345 


C 

C 

C 

C 


I 

I 


( 

I 

i 



The following computer programs and data files are recorded on this 
magtape. 


The order 
in magtape 

Descri ption 

Name suggested 
in IBM system 

1 

Particle trajectory computer program for 
axi symmetric inlet case. 

SOURCE. TRAJ3 

.2 

. 

Axisymmetric potential flow computer program (EOD) 

SOURCE. EOD 

3 

Axi symmetric COMBYN computer program 

SOURCE. COMBYN 

4 

Input data file (Tape 05) for geometry generation 
program SCIRCL (Axisymmetric inlet) 

DAT.AI005S 

5 

Input data file Tape 17 (renamed on Tape 05) 
for EOD 

DAT.AI005Y 

6 

Input data file (Tape 05) for COMBYN. 
(Axisymmetric inlet with a = 0°, M=0.4) 

DAT.AI005C 

7 

Input data file (Tape 05) for COMBYN 
(Axisymmetric inlet with a = 30°, M = 0.4) 

DAT.AI305C 

8 

Input data file (Tape 02) for Trajectory program, 
SOURCE. TRA03. (Axisymmetric inlet 
with a = 0°, M = 0.4) 

DAT.AI002T 




















9 

Input data file (Tape 02) for Trajectory program, 
SOURCE. TRAJ3. (Axi symmetric inlet 
with a = 30°, M=0.4) 

DAT.AI302T 

10 

Two-Dimensional potential flow computer program, 
Y24. (Original program). 

1 

SOURCE. Y240 

11 

Two-dimensional COMBIN computer program. (Original 
program) 

SOURCE. COMBI NO 

12 

Two-Dimensional COMBIN computer program. (Original 
program). This program is almost the same as 1 1 . 
We didn't use this one. 

Up to you. 

13 

Axi symmetric potential flow computer program, EOD. 
(Original program). 

SOURCE. EODO 

14 

Axi symmetric COMBYN computer program. (Original 
;■ program) . 

j 

SOURCE. COMBYNO 
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